c c ***************************************************************** c * * c * ROPA * c * program for rate-of-production analysis * c * * c ***************************************************************** c c ROPA is a member of KINAL , a program package for the c analysis of kinetic reaction mechanisms c c the following subroutines belong to ROPA main : c c EQS ORDER PIJ FUN JFY RTIME c ----------------------------------------------------------------- c c mode 0 Contribution of reaction steps to the c rate of production of species. c The least still listed contribution is greater than c 1/1000 of the most significant contribution. c c 1 The rate-of-production information is provided in c an abbreviated form: c c - The numbers of producing reactions are before the c slash (/), the numbers of consuming reactions are c after the slash. c - Both the producing and consuming reactions are in c descending order of magnitude. c - Reactions, having contribution less than 1/3 but c greater than 1/10 than that of the greatest c consuming/producing reactions, are given in parentheses. c - Reactions, having contribution less than 1/10 than c that of the greatest consuming/producing reactions, c are omitted. c - All the consuming reactions are omitted, if the c contribution of the greatest consuming reaction is less c than 1/100 than that of the greatest producing reaction. c All the producing reactions are omitted, if the c contribution of the greatest producing reaction is less c than 1/100 than that of the greatest consuming reaction. c c These limits (rlim1= 3., rlim2= 10., rlim3= 100.) c can be redefined altering the first executable c statements of the program. c c 2 Lifetime of species. c c ny - number of species ( max 50 ) c nr - number of reactions ( max 90 ) c nt - number of time points ( max 60 ) c nto - number of "observed" time points c lis - if lis.eq.0 listing of reactions is forbidden c c int - identification of "observed" time points in increasing order c i refers to the ith data block or to the (i+1)th time point c c formula - identification of species c y - concentrations c p - parameters c ta - time points ( ta(1) is the initial time ) c c c Warning: This program (mode 1) uses the character concatenation c operator (//). Some Fortran compilers may not know c this statement. c c ----------------------------------------------------------------- c c character*3 cnum(nr) c character*8 formula(ny) c dimension y(ny),dy(ny),ta(nt),int(nt),ip(nr),b(nr) c dimension fy(ny,ny),Pj(nr) 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*1 c, cj1, cj2 character*3 list, slash, lparen, rparen, cnum(90) character*8 formula(50),fal,comp character*10 datum,times character*30 fin,fout character*32 chm(3) character*80 rem dimension y(50),dy(50),ta(60),int(60),ip(90),b(90) dimension ifv1(20),ifv2(20),itv1(20),itv2(20),fy(50,50),Pj(90) 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 rlim1= 3. rlim2= 10. rlim3= 100. c comp= 'data No.' chm(1)= ' List of contributions ' chm(2)= ' The most significant reactions ' chm(3)= ' Lifetimes of species ' do 1 j=1,90 j1= j/10 j2= j-10*j1 cj1=char(j1+48) if (j1.eq.0) cj1=' ' cj2=char(j2+48) 1 cnum(j)= ' ' // cj1 // cj2 slash= ' / ' lparen= ' (' rparen= ') ' 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,nto,mode,lis read(5,100) (int(i),i=1,nto) 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) c list=' no' if (lis.ne.0) list='yes' c if (mode.lt.0.or.mode.gt.2) then write(6,229) mode stop endif c write(6,200) rem,chm(mode+1) c write(6,201) ny,nr,nt,nto,list c if (nto.gt.nt) then write(6,220) stop endif c do 2 i=1,nto if(int(i).lt.1.or.int(i).gt.nt) then write(6,202) int(i) stop endif if (i.eq.1) goto 2 if (int(i).le.int(i-1)) then write(6,222) stop endif 2 continue c write(6,204) (i,p(i),i=1,nr) write(6,205) (i,ta(i),i=1,nt) c nisa= 0 nisb= 0 3 read(5,104) k,l,ns if(k.eq.0) goto 4 if(k.gt.nr.or.l.gt.ny.or.ns.gt.3) write(6,207) k,l,ns if(k.le.0.or.l.le.0) write(6,207) k,l,ns if(ns.eq.0) ns = 1 nisa= nisa + 1 if (nisa.gt.nca) then write(6,227) nca stop endif isa(1,nisa) = k isa(2,nisa) = l isa(3,nisa) = ns goto 3 4 read(5,105) k,l,sn if(k.eq.0) goto 5 if(k.gt.nr.or.l.gt.ny.or.sn.gt.3.) write(6,208) k,l,sn if(k.le.0.or.l.le.0) write(6,208) k,l,sn nisb = nisb + 1 if (nisb.gt.ncb) then write(6,228) ncb stop endif isb(1,nisb) = k isb(2,nisb) = l sb(nisb) = sn goto 4 5 continue c if (lis.ne.0.and.nr.gt.7) write(6,209) if (lis.ne.0) call eqs(formula) c call rtime(datum,times) write(6,210) datum,times write(6,209) c------------------------------------------------------ do 51 it=1,nto 6 read (5,101) fal if (fal.eq.comp) goto 7 goto 6 7 read (5,106) ic if(ic.eq.int(it)) goto 8 goto 6 c 8 continue read (5,107) t,(y(i),i=1,ny) c write(6,203) int(it),t write(6,217) (formula(i),y(i),i=1,ny) c if (mode-1) 60, 70, 80 60 continue c c list of contributions c call fun(y,dy) do 61 ii=1,ny call pij(ii,y,Pj) call order(nr,ip,Pj,b) write(6,206) formula(ii), dy(ii) do 62 j=1,nr if( dabs(b(j))*1000. .lt. dabs(b(1)) ) goto 61 write(6,211) ip(j),b(j) 62 continue 61 continue goto 51 c 70 continue c c most significant reactions of a species c write(6,216) do 50 ii=1,ny call pij(ii,y,Pj) call order(nr,ip,Pj,b) c irp= 0 irm= 0 if1= 0 if2= 0 it1= 0 it2= 0 if (b(1)) 10, 53, 20 10 irm= 1 do 11 j=2,ny 11 if (b(j).gt.0.) goto 12 goto 30 12 fmax= b(1) tmax= b(j) if (tmax*rlim3.gt.-fmax) irp=1 goto 30 c 20 irp=1 do 21 j=2,ny 21 if (b(j).lt.0.) goto 22 goto 30 22 tmax= b(1) fmax= b(j) if (-fmax*rlim3.gt.tmax) irm=1 c c 30 do 52 j=1,20 if (b(j)) 31, 52, 40 31 if (b(j)*rlim1.lt.fmax) then if1= if1+1 ifv1(if1)= ip(j) elseif (b(j)*rlim2.lt.fmax) then if2= if2+1 ifv2(if2)= ip(j) endif goto 52 40 if (b(j)*rlim1.gt.tmax) then it1= it1+1 itv1(it1)= ip(j) elseif (b(j)*rlim2.gt.tmax) then it2= it2+1 itv2(it2)= ip(j) endif 52 continue c if (irm.eq.0) then if1=0 if2=0 endif if (irp.eq.0) then it1=0 it2=0 endif c 53 write(6,212) formula(ii),' ',(cnum(itv1(j)),j=1,it1),lparen, 1 (cnum(itv2(j)),j=1,it2),rparen,slash, 2 (cnum(ifv1(j)),j=1,if1),lparen, 3 (cnum(ifv2(j)),j=1,if2),rparen c 50 continue goto 51 c 80 continue c c lifetime of species c write(6,214) call jfy(ndy,y,fy) do 81 i=1,ny if (dabs(fy(i,i)).lt. 1.d-50) then Pj(i)= 9.999d+49 else Pj(i)= -1./fy(i,i) endif 81 continue call order(ny,ip,Pj,b) do 82 i=1,ny 82 write(6,215) i,formula(ip(i)),b(i) 51 continue call rtime(datum,times) write(6,210) datum,times write(6,209) 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) c 200 format(1x,a80//1x,a32) 201 format(//,3x,'number of species : ',i3/3x, 1 'number of reactions : ',i3/3x, 2 'number of time points : ',i3//3x, 3 'No of obs.time points : ',i3//3x, 4 'list of reactions : ',a3/) 202 format (/2x,'Error : unrecognisable time point reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 203 format(//' Time No: ',i3,' Time : ',1pe12.3/) 204 format (/2x,' Rate coefficients :'//100(4(2x,i5,1h.,1pe12.3)/)) 205 format (/2x,' Time points :'//100(4(2x,i5,1h.,1pe12.3)/)) 206 format (//1x,a8,' Rate : ',1pe12.3,' No Contribution:'/) 207 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i 3,' species :',i3,' stoich.coeff. :',i3) 208 format (/2x,' Warning ! Strange input data :'/ 2 2x,' reaction :',i3,' species :',i3,' stoich.coeff. :',f5.2) 209 format(1h1) 210 format(//20x,' Date : ',a10,' Time : ',a10//) 211 format(35x,i3,1pe15.3) 212 format(2x,a8,23a3) 214 format(/2x,'Lifetime of species : '//) 215 format(2x,i3,'. ',a8,1pe15.3) 216 format(//' Species: Significant reactions:'/) 217 format (/2x,' Concentrations :'// 100(3(2x,a8,2h =,1pd13.5)/)) 220 format (/2x,'Error : nto greater than nt'// 1 2x,' --- PROGRAM TERMINATED ---'///) 222 format (/2x,'Error : time point references are not '/ 1 /2x,' in increasing order '/// 2 /2x,' --- PROGRAM TERMINATED ---'///) 227 format(//2x,'Error : redimension isa of st common in main '// 1 2x,' and in the following subroutines '// 2 2x,' EQS , FUN , JFY , JFP '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '/ //) 228 format(//2x,'Error : redimension sb,isb of st common in main'// 1 2x,' and in the following subroutines '// 2 2x,' EQS , FUN , JFY , JFP '// 3 2x,' new dimension must be greater than',i4 // 4 2x,' --- PROGRAM TERMINATED --- '///) 229 format(//2x,'Error : ',i5,' is not existing mode '// 1 2x,' --- PROGRAM TERMINATED --- '///) c close (6, status = 'keep' ) stop 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 * ORDER puts into order the elements of a vector * c * according to their absolute value * c * * c **************************************************************** c c a - IN: vector of n elements (remains unchanged) c b - OUT: ordered vector c n - IN: dimension of vectors c ip - OUT: integer vector of n elements a(ip(i)) = b(i) c e - working area c c ---------------------------------------------------------------- c c dimension e(n) c subroutine order (n,ip,a,b) implicit integer*2(i-n) implicit real*8(a-h,o-z) dimension ip(1),a(1),b(1),e(100) do 1 i=1,n 1 e(i) = dsign(1.d+00,a(i)) c c rh must be greater then the greatest element of a rh = 1.d+50 do 2 i=1,n c c rm must be smaller then the smallest element of a rm = 1.d-50 do 3 j=1,n if(dabs(a(j)).ge.rm.and.dabs(a(j)).le.rh) goto 4 goto 3 4 if(i.eq.1) goto 5 do 6 k=2,i if(j.eq.ip(k-1)) goto 3 6 continue 5 rm=dabs(a(j)) im=j 3 continue b(i)=rm * e(im) rh=rm ip(i) = im 2 continue return end c c *** KINAL ****************************************************** c * * c * Pij matrix computation * c * * c **************************************************************** c c isp - IN: i c y - IN: actual concentrations c Pj - OUT: i-th row of Pij matrix c An element of this matrix contains the contribution c of the j-th reaction to the rate of production of c i-th species. 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 pij(isp,y,Pj) 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 y(1),Pj(1),r(90) c do 1 k=1,nr Pj(k)= 0. 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 is=1,nisa if(isp.eq.isa(2,is)) 1 Pj(isa(1,is))= - dble(isa(3,is)) * r(isa(1,is)) 3 continue do 4 is=1,nisb if(isp.eq.isb(2,is)) 1 Pj(isb(1,is))= Pj(isb(1,is)) + sb(is) * r(isb(1,is)) 4 continue 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 * 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