c c ***************************************************************** c * * c * YRED * c * program for variable elimination * c * * c ***************************************************************** c c YRED is a member of KINAL , a program package for the c analysis of kinetic reaction mechanisms c c The following subroutines belong to YRED main program: c c EQS ORDER FUN JFY RTIME c ----------------------------------------------------------------- c c ny - number of species ( max 50 ) | in the c nr - number of reactions ( max 90 ) | original c nt - number of time points ( max 60 ) | data file c c nyo - number of important species c c lis - if lis.eq.0 listing of reactions is forbidden c c iny - serial numbers of important species c int - serial number of "observed" 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 ndy - leading dimension of matrix fy c c ----------------------------------------------------------------- c c character*8 formula(ny) c dimension ta(nt),y(ny),dy(ny),iny(ny),fy(ny,ny) c dimension ip(ny),b(ny),bb(ny),ifn(ny),ifns(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*1 c,cm(2) character*3 list character*8 formula(50),fal,comp character*10 datum,times character*14 fin,fout character*25 chm(3) character*80 rem dimension ta(60),y(50),dy(50),iny(50),fy(50,50) dimension ip(50),b(50),bb(50),ifn(50),ifns(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.' cm(1)=' ' cm(2)='*' 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,nyo,lis read(5,100) (iny(i),i=1,nyo) read(5,101) (formula(i),y(i),i=1,ny) read(5,102) (p(i),i=1,nr) read(5,102) h read(5,102) (ta(i),i=1,nt) c list=' no' if (lis.ne.0) list='yes' c write(6,200) rem write(6,201) ny,nr,nt,nyo,list c if (nyo.gt.ny) then write(6,218) stop endif c do 2 i=1,nyo if(iny(i).lt.1.or.iny(i).gt.ny) then write(6,203) iny(i) stop endif 2 continue write(6,206) (formula(iny(i)),i=1,nyo) write(6,204) (i,p(i),i=1,nr) c write(*,223) read(*,*) int if (int.lt.1.or.int.gt.(nt-1)) then write(6,202) stop endif 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 c 6 read (5,101) fal if (fal.eq.comp) goto 7 goto 6 7 read (5,106) ic if(ic.eq.int) goto 8 goto 6 c 8 continue read (5,107) t,(y(i),i=1,ny) c call fun(y,dy) write (6,211) int,t write(6,217)(formula(i), y(i),i=1,ny) write(6,230)(formula(i),dy(i),i=1,ny) c do 9 i=1,ny 9 ifn(i)= 0 do 11 i=1,nyo 11 ifn(iny(i))= 1 call jfy(ndy,y,fy) icy= 0 c c begin of iteration loop c 30 continue icy= icy + 1 write(6,231) icy write(6,221) ( formula(iny(i)),i=1,nyo) c do 10 j=1,ny b(j)=0. do 10 k=1,nyo l=iny(k) if(dabs(dy(l)).lt.1.d-100) goto 10 b(j)= b(j)+(fy(l,j)*y(j)/dy(l))**2 10 continue c call order(ny,ip,b,bb) write(6,232) do 14 i=1,ny if (bb(i).lt.1.d-30) goto 14 bb(i)= dlog10(bb(i)) write(6,233) i,formula(ip(i)),cm(ifn(ip(i))+1),bb(i) 14 continue c 20 write(*,224) read (*,*) cl write(6,234) cl do 16 i=1,ny 16 if(b(i).gt.10.**cl) ifns(i)=1 c c list of new necessary species c in= 0 do 17 i=1,ny if(ifn(i).eq.0.and.ifns(i).eq.1) then in= in + 1 iny(in) = i endif 17 continue if (in.eq.0) goto 23 c write(6,235) ( formula(iny(i)),i=1,in) do 24 i=1,ny ifn(i)= ifn(i)+ifns(i) 24 if (ifn(i).ne.0) ifn(i)=1 c c list of necessary species c nyo= 0 do 25 i=1,ny if(ifn(i).gt.0) then nyo= nyo + 1 iny(nyo) = i endif 25 continue c goto 30 c 23 write(6,236) write(*,236) 12 write(*,237) read(*,109) c if(c.eq.'y') goto 22 if(c.eq.'n') goto 20 goto 12 c 22 continue c c final inventory c write(6,222) write(6,221) ( formula(iny(i)),i=1,nyo) nyo= 0 do 26 i=1,ny if (ifn(i).eq.0) then nyo= nyo + 1 iny(nyo) = i endif 26 continue write(6,225) ( formula(iny(i)),i=1,nyo) if (nyo.eq.0) write(6,226) c 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) 109 format (1a1) c 200 format(1x,a80) 201 format(//3x, 'number of species : ',i3/3x, 1 'number of reactions : ',i3/3x, 2 'number of time points : ',i3///3x, 3 'number of important species : ',i3//3x, 4 'list of reactions : ',a3/) 202 format (/2x,'Error : unrecognisable time point reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 203 format (/2x,'Error : unrecognisable obs. conc. reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 204 format (/2x,' Rate coefficients :'//100(4(2x,i5,1h.,1pe12.3)/)) 206 format (//2x,' Important species :',4(2x,a8)//10(7(2x,a8)//)) 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 (/' Reading data No.',i3, 1 ' Time : ',1pe15.3) 217 format (/2x,' Concentrations :'// 100(3(2x,a8,2h =,1pd13.5)/)) 218 format (/2x,'Error : nyo greater than ny'// 1 2x,' --- PROGRAM TERMINATED ---'///) 221 format (//1x,'Necessary species :',5(2x,a8)//10(7(2x,a8)//)) 222 format (/////2x,' **************** ITERATION TERMINATED' 1 ' ****************'//) 223 format(//' No. of time point ?') 224 format(//' Threshold value ?') 225 format (//1x,'Redundant species :',5(2x,a8)//10(7(2x,a8)//)) 226 format ( 1x,' none') 227 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 --- '/ //) 228 format(//2x,'Error : redimension sb,isb 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 --- '///) 230 format (/2x,' Production rates of species :' 1 // 100(3(2x,a8,2h =,1pd13.5)/)) 231 format (/////3x,'**************** Cycle',i3, 1 ' ****************') 232 format(/3x,'Estimated effect of species for the rate of' 1 ' necessary species (log units) :'//) 233 format(2x,i3,3x,a8,2x,a1,3x,f6.2) 234 format(////2x,'Threshold value :',2x,f6.2) 235 format (//2x,' New necessary species :'// 10(7(2x,a8)/)) 236 format (//2x,' No new necessary species !') 237 format (//2x,' Will you stop the iteration ? y/n ') 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 * 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