c----------------------------------------------------------------------- c Example program for GEARBI (revised version). c ODE system from diurnal kinetics transport PDE in 2 dimensions. c c This is the version of 2 August 2002. c This version is in double precision. c c The problem comes from a pair of coupled kinetics-diffusion c equations simulating the stratospheric chemistry of oxygen and c ozone in two space dimensions. The two species, O and O3, c have concentrations in moles/cc and undergo diurnal kinetics and c diffusion in the x and z directions, with diffusion coefficients c that are constant for the x diffusion and variable for z. c The region is 0 .le. x .le. 20, 30 .le. z .le. 50. c Zero flux boundary conditions are posed on all boundaries. c The equations are discretized with central differences on a uniform c mx by mz mesh to give the ODE system. The variables are ordered c first by species, then by x coordinate, then by z coordinate. c Thus the size parameters are mp = 2, mq = mx*mz, and n = 2*mx*mz. c The model time period is one day (in sec). c----------------------------------------------------------------------- implicit double precision (a-h,o-z) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 common /gear9/ hused, nqused, nstep, nfe, nje, nii dimension y0(800), fly(2), matp(2) data lout/6/ c mx and mz should be set to even numbers here. data mx/10/, mz/10/, dch/4.0d-6/, 1 q1/1.63d-16/, q2/4.66d-16/, a3/22.62d0/, a4/7.601d0/, 2 c3/3.7d16/, pi/3.1415926535898d0/, halfda/4.32d4/ c open (unit=6, file='kout', status='new') c npts = mx*mz mx2 = mx*2 midy = mz*mx - mx dx = 20.0d0/(mx-1) dz = 20.0d0/(mz-1) cox = dch/dx**2 coz = 1.0d0/dz**2 freq = pi/halfda nout = 12 tlast = 2.0d0*halfda floor = 1.0d2 c n = 2*npts t0 = 0.0d0 h0 = 2.0d-4 call yinit (y0) tout = 7200.0d0 eps = 1.0d-5 fly(1) = floor fly(2) = floor mf = 21 indx = 1 matp(1) = 2 matp(2) = npts c write(lout,20)mx,mz,n,t0,tlast,floor,eps,mf,h0 20 format('Example for GEARBI: 2-D diurnal kinetics-transport', 1 //' 2 species.. c1 = O and c2 = O3', 2 //' mesh is mx =',i3,' by mz =',i3,10x,' n =',i4, 3 //' t limits =',2e15.7, 4 //' semi-relative error control, floor =',e10.2, 5 //' eps =',e10.1,' mf =',i3,' h0 =',e12.3///) c do 80 iout = 1,nout call drivbi (n, t0, h0, y0, tout, eps, fly, mf, indx, matp) write(lout,70)tout,nstep,nfe,nje,nqused,h0 write(lout,71)y0(1),y0(2),y0(midy-1),y0(midy),y0(n-1),y0(n) 70 format('At t ='e12.4,' nstep =',i5,' nfe =',i5,' nje =',i4, 1 ' nq =',i2,' h =',e12.4) 71 format(20x,'At bottom left: c1, c2 =',2e12.4/ 1 20x,'At bottom right: c1, c2 =',2e12.4/ 2 20x,'At middle: c1, c2 =',2e12.4) if (indx .eq. 0) go to 78 write(lout,76)tout 76 format(//' final time reached =',e12.4//) stop 78 tout = tout + 7200.0d0 80 continue write(lout,85)nstep,nfe,nje,nii 85 format(//' final statistics..',i7,' steps,', 1 i8,' f-s,',i6,' j-s,',i7,' inner iters.') c mp = 2 na = 2 lwork = 6*n + 5*n + mp*mp*npts + 2*na*npts + mp write(lout,86)lwork 86 format(' total work space =',i8) c stop c end of main program for example problem. subroutines yinit, diffun, c pdbd, aset, offdi, diurn, chemr, and function routine dcv follow. end subroutine yinit (y0) c----------------------------------------------------------------------- c this routine computes and loads the vector of initial values. c the initial distribution for both species is a mildly peaked c function in x and z, consistent with the boundary conditions. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension y0(*) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c do 20 k = 1,mz z = 30.0d0 + (k-1)*dz z1 = 0.1d0*(z - 40.0d0) z1 = z1**2 cz = 1.0d0 - z1 + 0.5d0*z1**2 iz = 2*mx*(k-1) do 10 j = 1,mx x = (j-1)*dx x1 = 0.1d0*(x - 10.0d0) x1 = x1**2 cx = 1.0d0 - x1 + 0.5d0*x1**2 iy = iz + 2*(j-1) cxcz = cx*cz y0(iy+1) = 1.0d6*cxcz y0(iy+2) = 1.0d12*cxcz 10 continue 20 continue return c------------ end of subroutine yinit -------------------------------- end subroutine diffun (neq, t, y, ydot) c----------------------------------------------------------------------- c this routine computes the vector of right-hand sides of the ode c system, and returns it in ydot. diffun calls diurn and chemr for c the diurnal kinetics rates. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension y(neq), ydot(neq) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c call diurn(t) c do 90 k = 1,mz zl = 30.0d0 + (k - 1.5d0)*dz zu = zl + dz czdl = coz*dcv(zl) czdu = coz*dcv(zu) iz = mx2*(k-1) do 80 j = 1,mx iy1 = iz + 2*(j-1) + 1 iy2 = iy1 + 1 c1 = y(iy1) c2 = y(iy2) call chemr(c1,c2,rkin1,rkin2) if (k .ne. 1) go to 10 dcvu1 = y(iy1+mx2) - c1 dcvu2 = y(iy2+mx2) - c2 dcvl1 = -dcvu1 dcvl2 = -dcvu2 go to 30 10 continue if (k .ne. mz) go to 20 dcvl1 = c1 - y(iy1-mx2) dcvl2 = c2 - y(iy2-mx2) dcvu1 = -dcvl1 dcvu2 = -dcvl2 go to 30 20 continue dcvl1 = c1 - y(iy1-mx2) dcvl2 = c2 - y(iy2-mx2) dcvu1 = y(iy1+mx2) - c1 dcvu2 = y(iy2+mx2) - c2 30 continue if (j .ne. 1) go to 40 dchr1 = y(iy1+2) - c1 dchr2 = y(iy2+2) - c2 dchl1 = -dchr1 dchl2 = -dchr2 go to 60 40 continue if (j .ne. mx) go to 50 dchl1 = c1 - y(iy1-2) dchl2 = c2 - y(iy2-2) dchr1 = -dchl1 dchr2 = -dchl2 go to 60 50 continue dchl1 = c1 - y(iy1-2) dchl2 = c2 - y(iy2-2) dchr1 = y(iy1+2) - c1 dchr2 = y(iy2+2) - c2 60 continue ydot(iy1) = (czdu*dcvu1 - czdl*dcvl1) 1 + cox*(dchr1 - dchl1) + rkin1 ydot(iy2) = (czdu*dcvu2 - czdl*dcvl2) 1 + cox*(dchr2 - dchl2) + rkin2 80 continue 90 continue return c------------ end of subroutine diffun ------------------------------- end subroutine pdbd (t, yi, i, di, mp) c----------------------------------------------------------------------- c this routine computes the i-th 2 by 2 diagonal block of the jacobian c matrix and loads it into di. from the block index i, pdbd gets the c z coordinate index k in order to compute diffusion coefficients. c this routine assumes that subroutine diurn has just been called c at the same value of t. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension yi(2), di(2,2) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c k = (i-1)/mx + 1 zl = 30.0d0 + (k - 1.5d0)*dz zu = zl + dz czdl = coz*dcv(zl) czdu = coz*dcv(zu) tdiag = -(czdl + czdu + 2.0d0*cox) c1 = yi(1) c2 = yi(2) di(1,1) = -q1*c3 - q2*c2 + tdiag di(1,2) = -q2*c1 + q4 di(2,1) = q1*c3 - q2*c2 di(2,2) = -q2*c1 - q4 + tdiag return c------------ end of subroutine pdbd --------------------------------- end subroutine aset (t, y, aa, mq) c----------------------------------------------------------------------- c this routine computes the transport coefficents in the off-diagonal c blocks of the jacobian, and loads them into the mq by 4 array aa. c the off-diagonal block in block position (k,l) is a multiple a(k,l) c of the 2 by 2 identity matrix. the nonzero coefficients a(k,l) c are on the four lines k = l + 1, k = l - 1, k = l + mx, k = l - mx. c these are stored in columns 1 to 4 (respectively) of the array aa, c with a(k,l) stored in aa(k,m) for the appropriate value of m. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension y(*), aa(mq,4) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c mxm1 = mx - 1 cox2 = 2.0d0*cox iz = 0 do 90 k = 1,mz zl = 30.0d0 + (k - 1.5d0)*dz zu = zl + dz czdl = coz*dcv(zl) czdu = coz*dcv(zu) tzl = czdl if (k .eq. 1) tzl = 0.0d0 if (k .eq. mz) tzl = czdl + czdu tzu = czdu if (k .eq. 1) tzu = czdu + czdl if (k .eq. mz) tzu = 0.0d0 iq = iz + 1 aa(iq,1) = 0.0d0 aa(iq,2) = cox2 aa(iq,3) = tzl aa(iq,4) = tzu do 80 j = 2,mxm1 iq = iz + j aa(iq,1) = cox aa(iq,2) = cox aa(iq,3) = tzl aa(iq,4) = tzu 80 continue iq = iz + mx aa(iq,1) = cox2 aa(iq,2) = 0.0d0 aa(iq,3) = tzl aa(iq,4) = tzu iz = iq 90 continue return c------------ end of subroutine aset --------------------------------- end subroutine offdi (z, x, iq, con, aa, mp, mq) c----------------------------------------------------------------------- c this routine multiplies the vector x by the off-diagonal blocks of c block row iq of the jacobian (as stored in array aa), and adds c con times the result to the vector z (of length mp). c----------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension z(*), x(*), aa(mq,4) common /pcom/ dummy(12), mx, idummy(2) c aa1 = aa(iq,1) if (iq .lt. 2 .or. aa1 .eq. 0.0d0) go to 15 j1 = mp*(iq - 2) ca = con*aa1 do 10 ip = 1,mp 10 z(ip) = z(ip) + ca*x(j1+ip) 15 continue aa2 = aa(iq,2) if (iq .ge. mq .or. aa2 .eq. 0.0d0) go to 25 j2 = mp*iq ca = con*aa2 do 20 ip = 1,mp 20 z(ip) = z(ip) + ca*x(j2+ip) 25 continue aa3 = aa(iq,3) if (iq .le. mx .or. aa3 .eq. 0.0d0) go to 35 j3 = mp*(iq - mx - 1) ca = con*aa3 do 30 ip = 1,mp 30 z(ip) = z(ip) + ca*x(j3+ip) 35 continue aa4 = aa(iq,4) if (iq+mx .gt. mq .or. aa4 .eq. 0.0d0) go to 45 j4 = mp*(iq + mx - 1) ca = con*aa4 do 40 ip = 1,mp 40 z(ip) = z(ip) + ca*x(j4+ip) 45 continue return c------------ end of subroutine offdi -------------------------------- end subroutine diurn (t) c----------------------------------------------------------------------- c this routine computes the 2 diurnal rate functions at a given t. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c s = dsin(freq*t) if (s .le. 0.0d0) go to 20 q3 = dexp(-a3/s) q4 = dexp(-a4/s) return 20 q3 = 0.0d0 q4 = 0.0d0 return c------------ end of subroutine diurn -------------------------------- end subroutine chemr (c1, c2, r1, r2) c----------------------------------------------------------------------- c this routine computes the 2 kinetics rates at one spatial point. c----------------------------------------------------------------------- implicit double precision (a-h,o-z) common /pcom/ q1,q2,q3,q4,a3,a4,freq,c3,dx,dz,cox,coz,mx,mz,mx2 c qq1 = q1*c1*c3 qq2 = q2*c1*c2 qq3 = q3*c3 qq4 = q4*c2 r1 = -qq1 - qq2 + 2.0d0*qq3 + qq4 r2 = qq1 - qq2 - qq4 return c------------ end of subroutine chemr -------------------------------- end double precision function dcv (z) c----------------------------------------------------------------------- c this routine computes the vertical diffusion coefficient at a given z. c----------------------------------------------------------------------- double precision z, arg arg = 0.2d0*z dcv = 1.0d-8*dexp(arg) return c------------ end of function routine dcv ---------------------------- end Example for GEARBI: 2-D diurnal kinetics-transport 2 species.. c1 = O and c2 = O3 mesh is mx = 10 by mz = 10 n = 200 t limits = 0.0000000E+00 0.8640000E+05 semi-relative error control, floor = 0.10E+03 eps = 0.1E-04 mf = 21 h0 = 0.200E-03 At t = 0.7200E+04 nstep = 200 nfe = 300 nje = 24 nq = 5 h = 0.2101E+03 At bottom left: c1, c2 = 0.1040E+05 0.2509E+12 At bottom right: c1, c2 = 0.4029E+05 0.9725E+12 At middle: c1, c2 = 0.1121E+05 0.2705E+12 At t = 0.1440E+05 nstep = 226 nfe = 335 nje = 25 nq = 5 h = 0.3960E+03 At bottom left: c1, c2 = 0.6503E+07 0.2521E+12 At bottom right: c1, c2 = 0.2485E+08 0.9695E+12 At middle: c1, c2 = 0.7327E+07 0.2843E+12 At t = 0.2160E+05 nstep = 244 nfe = 365 nje = 27 nq = 4 h = 0.4167E+03 At bottom left: c1, c2 = 0.2604E+08 0.2919E+12 At bottom right: c1, c2 = 0.8512E+08 0.1005E+13 At middle: c1, c2 = 0.2947E+08 0.3334E+12 At t = 0.2880E+05 nstep = 262 nfe = 397 nje = 28 nq = 5 h = 0.3475E+03 At bottom left: c1, c2 = 0.8542E+07 0.3318E+12 At bottom right: c1, c2 = 0.2665E+08 0.1040E+13 At middle: c1, c2 = 0.9790E+07 0.3806E+12 At t = 0.3600E+05 nstep = 350 nfe = 548 nje = 47 nq = 4 h = 0.7506E+02 At bottom left: c1, c2 = 0.1380E+05 0.3330E+12 At bottom right: c1, c2 = 0.4297E+05 0.1037E+13 At middle: c1, c2 = 0.1608E+05 0.3879E+12 At t = 0.4320E+05 nstep = 444 nfe = 686 nje = 61 nq = 3 h = 0.5492E+03 At bottom left: c1, c2 = -0.1477E-09 0.3339E+12 At bottom right: c1, c2 = -0.4587E-09 0.1033E+13 At middle: c1, c2 = -0.1732E-09 0.3942E+12 At t = 0.5040E+05 nstep = 446 nfe = 688 nje = 62 nq = 3 h = 0.5492E+04 At bottom left: c1, c2 = -0.1476E-07 0.3348E+12 At bottom right: c1, c2 = -0.4577E-07 0.1030E+13 At middle: c1, c2 = -0.1734E-07 0.3998E+12 At t = 0.5760E+05 nstep = 447 nfe = 689 nje = 62 nq = 3 h = 0.5492E+04 At bottom left: c1, c2 = 0.2570E-11 0.3358E+12 At bottom right: c1, c2 = 0.9791E-11 0.1026E+13 At middle: c1, c2 = 0.2370E-11 0.4049E+12 At t = 0.6480E+05 nstep = 448 nfe = 690 nje = 62 nq = 3 h = 0.5492E+04 At bottom left: c1, c2 = -0.5181E-13 0.3367E+12 At bottom right: c1, c2 = -0.2252E-12 0.1022E+13 At middle: c1, c2 = -0.3317E-13 0.4095E+12 At t = 0.7200E+05 nstep = 450 nfe = 692 nje = 62 nq = 4 h = 0.7456E+04 At bottom left: c1, c2 = -0.2241E-13 0.3376E+12 At bottom right: c1, c2 = -0.3307E-13 0.1019E+13 At middle: c1, c2 = -0.3947E-13 0.4138E+12 At t = 0.7920E+05 nstep = 451 nfe = 693 nje = 62 nq = 4 h = 0.7456E+04 At bottom left: c1, c2 = -0.1182E-11 0.3386E+12 At bottom right: c1, c2 = -0.3773E-11 0.1015E+13 At middle: c1, c2 = -0.1367E-11 0.4178E+12 At t = 0.8640E+05 nstep = 452 nfe = 716 nje = 66 nq = 1 h = 0.6911E+03 At bottom left: c1, c2 = 0.1112E-11 0.3395E+12 At bottom right: c1, c2 = 0.4007E-11 0.1011E+13 At middle: c1, c2 = 0.1109E-11 0.4216E+12 final statistics.. 452 steps, 716 f-s, 66 j-s, 822 inner iters. total work space = 3002