c c ***************************************************************** c * * c * DIFF * c * program for the integration of * c * stiff kinetic differential equations * c * & check of reduced models * c * * c ***************************************************************** c c DIFF is a member of KINAL , a program package for the c analysis of kinetic reaction mechanisms c c The following subroutines belong to DIFF main program: c c ROW4A FUN JFY SOL DEC EQS RTIME c c ----------------------------------------------------------------- c c ny - number of species ( max 50 ) c nr - number of reactions ( max 90 ) c nt - number of time points ( max 60 ) c ipasto - max number of steps between two printouts c formula - identification of species c h - initial stepsize c hmax - maximal stepsize c tol - error tolerance in each step c y - concentrations at time t c p - parameters c ta - time points ( ta(1) is the initial time ) c lis - if lis.eq.0 listing of reactions is forbidden c c mode 0 output onto disc is forbidden c 1 output onto disc c 2 check of a reduced model c c If mode=2 is used, the file 'red.dat' has to be present in c the current directory. It contains a remark for the reduced c model (a80) and then two numbers (2I5) in each line: c The number of reactions in ascending order and 1 or 0, that c shows if the reaction is present or not in the reduced model. c c ----------------------------------------------------------------- c c character*8 formula(ny) c dimension ta(nt),y(ny),fy(ny,ny),ys(nt,ny),del(ny) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c c ----------------------------------------------------------------- c c written by Tamas Turanyi c Department of Physical Chemistry c E”tv”s University c c Address: H-1518 Budapest-112, P.O.Box 32, Hungary c Phone : (36-1) 209-0555 ext. 1109 c Fax : (36-1) 209-0602 c E-mail : turanyi@garfield.chem.elte.hu c WWW : www-PhCh.chem.elte.hu c c ----------------------------------------------------------------- c c refer to: T. Turanyi, Comp.Chem., Vol. 14, pp. 253-254, 1990 c c ----------------------------------------------------------------- c implicit real*8(a-h,o-z) implicit integer*2 (i-n) character*3 list character*8 formula(50),comp,fal character*10 datum,times character*14 fin,fout character*30 chm(3) character*80 rem dimension ta(60),y(50),fy(50,50),ys(60,50),del(50) common /st / nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c ndy= 50 nca= 200 ncb= 300 c comp = 'data No.' chm(1)= ' Output onto disc is forbidden' chm(2)= ' Output onto disc ' chm(3)= ' Check of a reduced model ' c write(*,*) 'Name of data file ?' read(*,'(a)') fin write(*,*) 'Name of output file ?' read(*,'(a)') fout c open(5,file=fin,status='old') open(6,file=fout,status='new') c read(5,103) rem read(5,100) ny,nr,nt,mode,lis read(5,101) (formula(i),y(i),i=1,ny) read(5,102) (p(i),i=1,nr) read(5,102) h,hmax,tol read(5,102) (ta(i),i=1,nt) write(6,200) rem c if (mode.lt.0.or.mode.gt.2) then write (6,203) mode stop endif if (mode.ne.2) goto 9 open(4,file='red.dat',status='old') ns= 0 read(4,103) rem do 8 i=1,nr read(4,100) n,nred if (nred.ne.0) nred=1 p(i)= p(i)* dble(nred) 8 ns=ns+nred write(6,214) ns write(6,200) rem c 9 list=' no' if(lis.ne.0) list='yes' write(6,201) ny,nr,nt,list write(6,204) mode,chm(mode+1) write(6,205) (i,p(i),i=1,nr) write(6,206) h,hmax,tol write(6,207) (i,ta(i),i=1,nt) write(6,208) (formula(i),y(i),i=1,ny) c nisa = 0 nisb = 0 3 read(5,104) k,l,ns if(k.eq.0) goto 2 if(k.gt.nr.or.l.gt.ny.or.ns.gt.3) write(6,209) k,l,ns if(k.le.0.or.l.le.0) write(6,209) k,l,ns if(ns.eq.0) ns = 1 nisa = nisa + 1 if (nisa.gt.nca) then write(6,220) nca stop endif isa(1,nisa) = k isa(2,nisa) = l isa(3,nisa) = ns goto 3 2 read(5,105) k,l,sn if(k.eq.0) goto 1 if(k.gt.nr.or.l.gt.ny.or.sn.gt.3.) write(6,210) k,l,sn if(k.le.0.or.l.le.0) write(6,210) k,l,sn nisb = nisb + 1 if (nisb.gt.ncb) then write(6,221) ncb stop endif isb(1,nisb) = k isb(2,nisb) = l sb(nisb) = sn goto 2 1 continue c if (mode.ne.2) goto 12 do 11 i=1,(nt-1) 4 read (5,101) fal if (fal.eq.comp) goto 5 goto 4 5 read (5,106) ic if(ic.eq.i) goto 6 goto 4 6 continue c read (5,107) t,(ys(i,j),j=1,ny) c write (6,212) i,t write(6,208)(formula(j),ys(i,j),j=1,ny) 11 continue c 12 if (nr.gt.7.and.lis.ne.0) write(6,219) if (lis.ne.0) call eqs(formula) c call rtime(datum,times) write(6,211) datum,times c ipa= 0 ire= 0 nt= nt - 1 do 10 ict=1,nt t=ta(ict) 7 tend=ta(ict+1) ipas= 500 call row4a(ny,t,y,fy,tend,tol,hmax,h,ipas,irep) ipa= ipa+iabs(ipas) ire= ire+irep if(ipas.lt.0) goto 7 write(6,213) t write(6,208) (formula(k),y(k),k=1,ny) write(6,215) ipa,ire if (mode.ne.1) goto 14 write(5,101) comp write(5,106) ict write(5,107) t,(y(i),i=1,ny) 14 continue c if (mode.ne.2) goto 10 do 15 i=1,ny if(dabs(ys(ict,i)).lt.1.d-50) then del(i)= 0. else del(i)=(y(i)-ys(ict,i))/ys(ict,i)*100. endif 15 continue write(6,202) (formula(k),del(k),k=1,ny) 10 continue c call rtime(datum,times) write(6,211) datum,times c 100 format(16i5) 101 format(a8,2x,f10.2) 102 format (8f10.2) 103 format(a80) 104 format(10i5) 105 format(2i5,f5.2) 106 format (i3) 107 format (6x,6e12.5) 108 format (i3,3x,6e12.5,10(/6x,6e12.5)) c 200 format(1x,a80) 201 format(//,3x,'number of species :',i3/3x, 2 'number of reactions :',i3/3x, 3 'number of time points :',i3/3x, 4 'listing of reactions :',a3/) 202 format (/2x,' Deviations :'// 100(3(2x,a8,2h =,1pg13.5)/)) 203 format (/2x,'Error : ',i5,' is not existing mode '// 1 2x,' --- PROGRAM TERMINATED ---'///) 204 format(1x,'Mode : ',i2,2x,a30) 205 format (/2x,' Rate coefficients :'//100(4(2x,i5,1h.,1pe12.3)/)) 206 format(//,2x,' h :',1pd11.4/2x, 2 ' hmax :',1pd11.4/2x, 3 ' tol :',1pd11.4/2x) 207 format (/2x,' Time points :'//100(4(2x,i5,1h.,1pe12.3)/)) 208 format (/2x,' Concentrations :'// 100(3(2x,a8,2h =,1pd13.5)/)) 209 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i3,' species :',i3,' stoich.coeff. :',i3) 210 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i3,' species :',i3,' stoich.coeff. :',f5.2) 211 format(//20x,' Date : ',a10,' Time : ',a10//) 212 format(//2x,'No. of data set :',i3,6x,'Time :',1pe10.2) 213 format(/2x,' Time :',1pd13.4,5x/) 214 format(/5x,i3,' step reduced model') 215 format(//2x,' number of steps : ',i4/ 1 2x,' number of backsteps : ',i4) 219 format(1h1) 220 format(//2x,'Error : redimension isa of ST common in main '// 1 2x,' and in the following subroutines : '// 2 2x,' EQS , FUN , JFY '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '///) 221 format(//2x,'Error : redimension isb,sb of ST common in main'// 1 2x,' and in the following subroutines : '// 2 2x,' EQS , FUN , JFY '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '///) c close (6, status = 'keep' ) stop end c c *** KINAL ******************************************************* c * * c * ROW4A integrates stiff differential equations * c * * c ***************************************************************** c c A reliable Rosenbrock integrator for stiff differential equations c B.A. Gottwald , G. Wanner : Computing 26(4),355-360(1981) c c the original program is slightly modified c c ----------------------------------------------------------------- c c n - IN: number of variables ( max 50 ) c t - IN: initial time c OUT: final time c y - IN: initial y values c OUT: the same at final time c e - working area of minimum n*n elements c tend IN: end time c tol IN: relative error tolerance per step c (suggested value: 1.d-2 - 1.d-4) c hmax IN: maximal step size c h IN: initial step size c OUT: step size used by ROW4A c ipas IN: maximum number of steps c OUT: if ipas>0 number of succesful steps c if ipas<0 ROW4A failed to solve the problem c (try to increase ipas or to decrease tol) c irep OUT: number of repeated steps c c ----------------------------------------------------------------- c c dimension ynew(n),yold(n),dy(n) c dimension ak1(n),ak2(n),ak3(n),ak4(n) c dimension ip(n) c subroutine row4a(n,t,y,e,tend,tol,hmax,h,ipas,irep) implicit integer*2(i-n) implicit real*8(a-h,o-z) logical*2 end dimension y(1),ynew(50),yold(50),dy(50),e(n,1), 1 ak1(50),ak2(50),ak3(50),ak4(50),ip(50) c c---------------------------------------------------------------------- ipasto=ipas c---------------------initiate variables------------------------------- a21= .438000000000000d+00 a31= .938948678483428d+00 a32= .730795420615381d-01 c21=-.194347441894707d+01 c31= .416957530989189d+00 c32= .132396782072923d+01 c41= .151951325778448d+01 c42= .135370815030093d+01 c43=-.854151495257539d+00 b1 = .729044879960308d+00 b2 = .541069773272405d-01 b3 = .281599362440017d+00 b4 = .250000000000000d+00 e1 =-.190858871999474d-01 e2 = .255608791716455d+00 e3 =-.863816280897592d-01 e4 = .250000000000000d+00 ipas=0 irep=0 if(tend.eq.t) goto 55 if(tend.lt.t) goto 77 uround=dmax1(1.73d-18,tol*1.d-08) if(h.eq.0.d0) h=1.d-02 h=dabs(h) hold=0.d+00 end=.false. facmx=2.d+00 facmn=.125d+00 c c--------------------compute new step---------------------------------- 11 continue h=dmin1(h,(tend-t)) if(ipas.ge.ipasto) goto 77 c--------------------compute Jacobian----------------------------------- call fun(y,dy) 12 call jfy(n,y,e) c--------------compute matrix of linear system-------------------------- 17 do 200 i=1,n do 200 j=1,n 200 e(i,j)=-h*e(i,j)*0.395d+00 do 201 i=1,n 201 e(i,i)=1.d+00+e(i,i) call dec(n,e,ip,ierr) if(ierr.ne.0) h=h/2 if(ierr.ne.0) goto 12 c-------------------the stages----------------------------------------- do 203 i=1,n 203 ak1(i)=h*dy(i) call sol(n,e,ip,ak1) do 210 i=1,n 210 ynew(i)=y(i)+a21*ak1(i) call fun(ynew,dy) do 211 i=1,n 211 ak2(i)=h*dy(i)+c21*ak1(i) call sol(n,e,ip,ak2) do 220 i=1,n 220 ynew(i)=y(i)+a31*ak1(i)+a32*ak2(i) call fun(ynew,dy) do 221 i=1,n 221 ak3(i)=h*dy(i)+c31*ak1(i)+c32*ak2(i) call sol(n,e,ip,ak3) do 231 i=1,n 231 ak4(i)=h*dy(i)+c41*ak1(i)+c42*ak2(i)+c43*ak3(i) call sol(n,e,ip,ak4) c------------------------new solution------------------------------ do 240 i=1,n 240 ynew(i)=y(i)+b1*ak1(i)+b2*ak2(i)+b3*ak3(i)+b4*ak4(i) tnew=t+h c--------------------compute error estimation----------------------- est=uround do 300 i=1,n s=e1*ak1(i)+e2*ak2(i)+e3*ak3(i)+e4*ak4(i) sk=dmax1(1.d-30,y(i),ynew(i)) 300 est=dmax1(est,dabs(s)/sk) c--------------------new step size------------------------------------ if(.not.end) *hnew=dmin1(hmax,h*dmin1(facmx,dmax1(facmn,.9e0*(tol/est)**.25e0))) c---------------------error test--------------------------------------- if(end) goto 51 if(est.gt.tol) goto 43 c-----------------------accept step------------------------------------- if(t.gt.tend) goto 48 facmn=.5d+00 if(facmx.ne..95d0) facmx=1.5d0 if(facmx.eq..95d0) facmx=0.9d0 ipas=ipas+1 do 41 i=1,n yold(i)=y(i) 41 y(i)=ynew(i) told=t t=tnew hold=h h=hnew write(*,100) ipas,t,h,est 100 format(5x,i5,' t= ',1pd10.3,' h= ',1pd10.3, 1 ' est= ',1pd10.3) if(t.eq.tend) return goto 11 c------------------test for necessary back-step-------------------------- 43 if(ipas.gt.0) irep=irep+1 if(hnew.ge.hold) goto 47 c--------------------------back step----------------------------------- h=hnew ipas=ipas-1 irep=irep+1 facmx=.95d0 44 t=told do 45 i=1,n 45 y(i)=yold(i) hold=0.d0 goto 11 c-----------------repeat last step only---------------------------------- 47 h=hnew facmx=1.d0 goto 11 c-----------------manage the final step------------------------------ 48 h=tend-told end=.true. goto 44 c----------------------successful return--------------------------- 51 t=tend do 53 i=1,n 53 y(i)=ynew(i) 55 continue h=hnew return c--------------------unsuccessful return----------------------------- 77 ipas=-ipas return end c c *** KINAL ****************************************************** c * * c * FUN computes the right hand side of kinetic diff. eq. * c * * c **************************************************************** c c y - IN: actual concentrations c dy - OUT: derivative of variables with respect to time c dy(i) = d y(i) / d t c r - working area c c ---------------------------------------------------------------- c c dimension r(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine fun (y,dy) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension y(1),dy(1),r(90) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c do 1 k=1,nr 1 r(k)= p(k) do 2 is=1,nisa 2 r(isa(1,is)) = r(isa(1,is))*y(isa(2,is))**isa(3,is) do 3 i=1,ny 3 dy(i) = 0. do 4 is=1,nisa 4 dy(isa(2,is))=dy(isa(2,is)) - dble(isa(3,is))*r(isa(1,is)) do 5 is=1,nisb 5 dy(isb(2,is))=dy(isb(2,is)) + sb(is)*r(isb(1,is)) return end c c *** KINAL ****************************************************** c * * c * JFY computes the Jacobian with respect to variables * c * * c **************************************************************** c c y - IN: actual concentrations c f - OUT: Jacobian matrix c fy(i,j) = d f(i) / d y(j) c ndy - IN: leading dimension of matrix fy c r - working area c c ---------------------------------------------------------------- c c dimension r(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine jfy(ndy,y,fy) implicit integer*2(i-n) implicit real*8(a-h,o-z) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) dimension fy(ndy,1),y(1),r(90) c do 1 i=1,ny do 1 j=1,ny 1 fy(i,j)=0. do 6 j=1,ny do 2 k=1,nr 2 r(k) = 0. do 3 is=1,nisa 3 if(isa(2,is).eq.j) r(isa(1,is)) = p(isa(1,is)) do 4 is=1,nisa if(r(isa(1,is)).eq.0.) goto 4 if(isa(2,is).eq.j) then r(isa(1,is)) = r(isa(1,is))*dble(isa(3,is)) 1 *y(isa(2,is))**(isa(3,is)-1) else r(isa(1,is)) = r(isa(1,is))*y(isa(2,is))**isa(3,is) endif 4 continue do 5 is=1,nisa if (r(isa(1,is)).eq.0.) goto 5 fy(isa(2,is),j)=fy(isa(2,is),j)-dble(isa(3,is))*r(isa(1,is)) 5 continue do 6 is=1,nisb if (r(isb(1,is)).eq.0.) goto 6 fy(isb(2,is),j)=fy(isb(2,is),j) + sb(is) * r(isb(1,is)) 6 continue return end c c *** KINAL ****************************************************** c * * c * SOL computes the solution of linear system A*X=B * c * * c **************************************************************** c c see C.B. Moler , Algorithm 423 , CACM 15 (1972) c subroutine sol(n,a,ip,b) integer*2 n,ip(1),nm1,k,kp1,m,i,kb,km1 double precision a(n,1),b(1),t if (n.eq.1) goto 50 nm1=n-1 do 20 k=1,nm1 kp1=k+1 m=ip(k) t=b(m) b(m)=b(k) b(k)=t do 10 i=kp1,n 10 b(i)=b(i)+a(i,k)*t 20 continue do 40 kb=1,nm1 km1=n-kb k=km1+1 b(k)=b(k)/a(k,k) t=-b(k) do 30 i=1,km1 30 b(i)=b(i)+a(i,k)*t 40 continue 50 b(1)=b(1)/a(1,1) return end c c *** KINAL ****************************************************** c * * c * DEC makes the Gaussian decomposition of matrix A * c * * c **************************************************************** c c see C.B. Moler , Algorithm 423 , CACM 15 (1972) c subroutine dec(n,a,ip,ier) integer*2 n,ip(1),ier,nm1,k,kp1,m,i,j double precision a(n,1),t ier=0 ip(n)=1 if(n.eq.1) goto 70 nm1=n-1 do 60 k=1,nm1 kp1=k+1 m=k do 10 i=kp1,n 10 if(dabs(a(i,k)).gt.dabs(a(m,k))) m=i ip(k)=m t=a(m,k) if(m.eq.k) goto 20 ip(n)=-ip(n) a(m,k)=a(k,k) a(k,k)=t 20 if (t.eq.0.d00) goto 80 t=1.d00/t do 30 i=kp1,n 30 a(i,k)=-a(i,k)*t do 50 j=kp1,n t=a(m,j) a(m,j)=a(k,j) a(k,j)=t if(t.eq.0.d00) goto 50 do 40 i=kp1,n 40 a(i,j)=a(i,j)+a(i,k)*t 50 continue 60 continue 70 k=n if (a(n,n).eq.0.d00) goto 80 return 80 ier=k ip(n)=0 return end c c *** KINAL ****************************************************** c * * c * EQS makes the list of reactions * c * * c **************************************************************** c c formula - IN: identification of species c subroutine eqs(formula) implicit integer*2 (i-n) implicit real*8 (a-h,o-z) character left(3) character*3 plus(10) character*8 blank,formula(1),aleft(3),aright(10) dimension right(10) c common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) data plus / 10*' + ' / c c the elements of isa ; isb,sb vectors are assumed c to be arranged in increasing order of reactions c iper = 6 blank =' ' jsa = 0 jsb = 0 do 1 k=1,nr il = 0 ir = 0 do 2 i=1,3 left(i) = ' ' 2 aleft(i)=blank do 3 i=1,10 right(i) = 0. 3 aright(i)=blank c 4 continue jsa= jsa + 1 if (jsa.gt.nisa) goto 7 if (isa(1,jsa).ne.k) goto 5 il= il + 1 if(il.gt.3) then il= 3 write(iper,100) k goto 7 endif aleft(il) = formula(isa(2,jsa)) if(isa(3,jsa).eq.2) left (il) = '2' if(isa(3,jsa).eq.3) left (il) = '3' goto 4 5 continue jsa = jsa - 1 7 continue jsb= jsb + 1 if (jsb.gt.nisb) goto 8 if (isb(1,jsb).ne.k) goto 6 ir= ir + 1 if(ir.gt.10) then ir= 10 write(iper,110) k goto 8 endif aright(ir) = formula(isb(2,jsb)) right(ir) = sb(jsb) goto 7 6 continue jsb = jsb - 1 8 continue c il1 = il - 1 ir3 = ir - 3 ir2 = min0(ir,2) ir1 = ir2 - 1 write(iper,120) k,(left(i),aleft(i),i=1,3), 1 (right(i),aright(i),i=1,ir2) write(iper,130) (plus(i),i=1,il1) write(iper,140) if(ir1.gt.0) write(iper,150) plus(1) if(ir .gt.2) then write(iper,160) (right(i),aright(i),i=3,ir) endif c 1 continue c 100 format(1x,'More then 3 species on the left side of eqs. No ',i3) 110 format(1x,'More then 10 species on the right side of eqs. No ',i3) 120 format(/1x,i3,'. ',3(a1,1x,a8,3x),2x,2(f4.2,1x,a8,3x)) 130 format(1h+,15x,3(a3,11x)) 140 format(1h+,41x,' --> ') 150 format(1h+,59x,a3,13x,a3) 160 format(2(/9x,4(' + ',f5.2,1x,a8))) 170 format(1h+,22x,8(a3,13x)) return end c c *** KINAL ****************************************************** c * * c * RTIME reads date and time * c * * c **************************************************************** c c This subroutine is machine and compiler dependent. c Recent version conforms to MS-FORTRAN 4.0 c subroutine rtime(cd,ct) implicit integer*2 (i-n) character*10 cd,ct c c date c call getdat(iy,im,id) iy1=iy/100 iy=iy-iy1*100 iy1=iy/10 iy2=iy-iy1*10 +48 iy1=iy1 +48 im1=im/10 im2=im-im1*10 +48 im1=im1 +48 id1=id/10 id2=id-id1*10 +48 id1=id1 +48 cd=char(im1)//char(im2)//'-'//char(id1)//char(id2)//'-' 1 //char(iy1)//char(iy2) c c time c call gettim(ih,im,is,ic) ih1=ih/10 ih2=ih-ih1*10 +48 ih1=ih1 +48 im1=im/10 im2=im-im1*10 +48 im1=im1 +48 is1=is/10 is2=is-is1*10 +48 is1=is1 +48 ct=char(ih1)//char(ih2)//':'//char(im1)//char(im2)//':' 1 //char(is1)//char(is2) return end