c c ***************************************************************** c * * c * SENS * c * program for computing the local sensitivity matrix * c * * c ***************************************************************** c c SENS is a member of KINAL , a program package for the c analysis of kinetic reaction mechanisms c c The following subroutines belong to SENS main program: c c EQS RMTXKI ROW4SP LINS FCNS SOLS AU c DEC SOL FUN JFY JFP JF2YP JF2YY c RTIME c c ----------------------------------------------------------------- c c ny - number of species ( max 50 ) c nr - number of reactions ( max 90 ) c np - number of parameters ( max 20 ) c nt - number of time points ( max 60 ) c inp - serial number of parameters respect to c sensitivity computation is required c 1 <= inp(i) <= nr refers to rate coefficients c (nr+1) <= inp(i) <= (nr+ny) refers to init. concentrations c ipasto - max number of steps between two printouts c ifo - number of last computed matrix c formula - identification of species c h - initial stepsize c hmax - maximal stepsize c c error tolerance in each step at the c tol - solution of diff. equations c tols - sensitivity computing c c y - concentrations and sensitivities 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 ----------------------------------------------------------------- c c n = ny * (np+1) c c character*8 formula(ny) c dimension ta(nt),y(n),inp(np),ink(nr),inc(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),fal,comp character*10 datum,times character*14 fin,fout character*80 rem dimension ta(60),y(1050),inp(20),ink(90),inc(50) common /st/ nisa,nisb,isa(3,200),isb(2,300),sb(300), 1 nr,ny,p(90) c nca= 200 ncb= 300 comp= 'data No.' 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,np,nt,ifo,lis read(5,100) (inp(i),i=1,np) read(5,101) (formula(i),y(i),i=1,ny) read(5,102) (p(i),i=1,nr) read(5,102) h,hmax,tol,tols 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,np,nt,ifo,list write(6,205) (i,p(i),i=1,nr) write(6,206) h,hmax,tol,tols write(6,207) (i,ta(i),i=1,nt) write(6,208) (formula(i),y(i),i=1,ny) c ik= 0 ic= 0 do 11 i=1,np if(inp(i).le.0) goto 12 if(inp(i).le.nr) goto 14 if(inp(i).le.(nr+ny)) goto 15 12 write(6,216) inp(i) stop 14 ik= ik + 1 ink(ik)= inp(i) goto 11 15 ic= ic + 1 inc(ic)= inp(i) - nr 11 continue write(6,217) if(ik.eq.0) write(6,203) if(ik.ne.0) write(6,218) (ink(i),i=1,ik) write(6,202) if(ic.eq.0) write(6,203) if(ic.ne.0) write(6,204) (formula(inc(i)),i=1,ic) 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 (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 ny1= ny+1 do 9 j=1,np do 9 i=1,ny k= j*ny+i y(k)= 0. if(inp(j).le.nr) goto 9 k= j*ny+inp(j)-nr y(k)= 1. 9 continue c if(ifo.le.0) goto 8 c 4 read (5,101) fal if (fal.eq.comp) goto 5 goto 4 5 read (5,106) ic if(ic.eq.ifo) goto 6 goto 4 6 continue c read (5,107) t,(y(i),i=1,ny) do 20 j=1,np read (5,107) (y(i),i=(j*ny+1),(j*ny+ny)) 20 continue c write (6,212) ifo,t write(6,208)(formula(i),y(i),i=1,ny) call rmtxki (ny,np,ny,y(ny1)) 8 continue c ipa= 0 ire= 0 ifo= ifo + 1 nt= nt - 1 do 10 ict=ifo,nt t=ta(ict) 7 tend=ta(ict+1) ipas= 500 call row4sp(np,ny,t,y,tend,tol,tols,hmax,h,ipas,irep,inp) 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 (np.gt.0) call rmtxki (ny,np,ny,y(ny1)) c write(5,101) comp write(5,106) ict write(5,107) t,(y(i),i=1,ny) do 21 j=1,np write(5,108) inp(j),(y(i),i=(j*ny+1),(j*ny+ny)) 21 continue 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 parameters :',i3/3x, 4 'number of time points :',i3/3x, 5 'first sens. matrix :',i3/3x, 6 'listing of reactions :',a3/) 202 format (/2x,' Sens. calculation for the initial', 2 ' conc. of following species :'/) 203 format (2x,' none') 204 format (2x,5(a8,2x)) 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, 4 ' tols :',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 (1h1///' The last sensitivity matrix No.',i3, 1 / ' time : ',1pe15.3) 213 format(/2x,' Time :',1pd13.4,5x/) 215 format(//2x,' number of steps : ',i4/ 1 2x,' number of backsteps : ',i4) 216 format (/2x,'Error : unrecognizable parameter reference: ',i3// 1 2x,' --- PROGRAM TERMINATED ---'///) 217 format (/2x,' Sens. calculation for the following' 2 ' rate coefficients :'/) 218 format (15(2x,i3)/) 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 , JFP , JF2YY , JF2YP '// 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 , JFP , JF2YY , JF2YP '// 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 * 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 * RMTXKI writes large real matrices * c * * c ***************************************************************** c c r - matrix to write c nd - leading dimension of r c n - number of rows of matrix r c m - number of columns of matrix r c subroutine rmtxki(n,m,nd,r) implicit real*8(a-h,o-z) implicit integer*2(i-n) dimension r(nd,m) nm=m/5+1 nl=m do 1 k=1,nm nk=min0 (nl,5) nn = (k-1)*5 + nk ni = (k-1)*5 + 1 write(6,10) (i,i=ni,nn) do 2 i=1,n 2 write (6,11) i,(r(i,j),j=ni,nn) nl=nl-5 if(nl.eq.0) return 1 continue 10 format(///,6x,10(11x,i3)/) 11 format(/2x,i3,2x,5(1pe14.4)) return end c c *** KINAL ******************************************************* c * * c * ROW4SP computes the local sensitivity coefficients * c * * c ***************************************************************** c c An extended ODE solver for sensitivity calculations c P. Valko and S. Vajda : Computers & Chemistry 8(4),255-271(1984) c c ROW4SP is a PC version of ROW4S c FCNS , LINS , SOLS and AU were also modified. c Due to modifications ROW4SP consumes considerable less memory space, c but is much slower then the original ROW4S. c Use the original version if you have enough memory space. c c ROW4SP is an extension of ROW4A ( see B.A. Gottwald and G. Wanner c Computing 26(4),355-360(1981) ) c c ----------------------------------------------------------------- c c np - IN: number of parameters, related to which c sensitivities are to be computed c at np=0 only the original ODE-s will c be solved - in that case JFP,JF2YY,JF2YP c can be dummy subroutines c ny - IN: number of dependent variables c t - IN: initial time c OUT: final time c y - IN: first ny elements are actual y values c next ny elements are sensitivities with c respect to the first parameter,etc. c OUT: the same at final time c tend - IN: end time c tol - IN: relative error tolerance for the ODE solution c (suggested value: 1.d-2 - 1.d-4) c tols - IN: relative error tolerance for sensitivities c (suggested value: 1.d0 - 1.d-2) c hmax - IN: maximal step size c h - IN: initial step size c OUT: step size suggested by ROW4SP c ipas - IN: maximal number of steps c OUT: if ipas>0 number of succesful steps c if ipas<0 ROW4SP failed to solve the problem c (try to increase ipas or to decrease tol or tols) c irep - OUT: number of back steps c inp - IN: integer vector containing the serial number c of parameters , related to which c sensitivities are to be computed c c ----------------------------------------------------------------- c c n = ny * (np+1) c nyy= ny * ny c c dimension dy(n),ynew(n),yold(n),ak1(n),ak2(n),ak3(n),ak4(n) c dimension e(nyy),fy(nyy),f2yp(nyy),ip(ny) c subroutine row4sp + (np,ny,t,y,tend,tol,tols,hmax,h,ipas,irep,inp) implicit real*8(a-h,o-z) implicit integer*2(i-n) logical*2 end dimension y(1),inp(1) common /w1/ dy(1050),ynew(1050),yold(1050), 1 ak1(1050),ak2(1050),ak3(1050),ak4(1050) common /w2/ e(2500),fy(2500),f2yp(2500),ip(50) c ------------------------------------------------------------- n= ny*(np+1) ipasto = ipas c------------------initialize 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) go to 55 if (tend.lt.t) go to 77 uround= dmax1(1.73d-18,tol*1.d-08) if (h.eq.0.d0) h=1.d-2 h= dabs(h) hold= 0.d0 end= .false. facmx= 2.d0 facmn= .125d0 c----------------------new step---------------------------- 11 continue h=dmin1(h,(tend-t)) if (ipas.ge.ipasto) goto 77 call fcns(np,ny,y,dy,fy,inp) call lins(np,ny,h,y,fy,e,f2yp,inp) call dec(ny,e,ip,ierr) if (ierr.ne.0) h=h/2.d0 if (ierr.ne.0) goto 11 c-----------------------the stages--------------------------------- do 203 i=1,n 203 ak1(i)= h*dy(i) call sols(np,ny,e,f2yp,ip,ak1) do 210 i=1,n 210 ynew(i)= y(i)+a21*ak1(i) call fcns(np,ny,ynew,dy,fy,inp) do 211 i=1,n 211 ak2(i) = h*dy(i)+c21*ak1(i) call sols(np,ny,e,f2yp,ip,ak2) do 220 i=1,n 220 ynew(i)= y(i)+a31*ak1(i)+a32*ak2(i) call fcns(np,ny,ynew,dy,fy,inp) do 221 i=1,n 221 ak3(i)= h*dy(i)+c31*ak1(i)+c32*ak2(i) call sols(np,ny,e,f2yp,ip,ak3) do 231 i=1,n 231 ak4(i)= h*dy(i)+c41*ak1(i)+c42*ak2(i)+c43*ak3(i) call sols(np,ny,e,f2yp,ip,ak4) c-------------------------new solution------------------------------- do 240 i=1,n 240 ynew(i)= y(i)+ 1 b1*ak1(i)+b2*ak2(i)+b3*ak3(i)+b4*ak4(i) tnew= t+h c----------------------error estimation------------------------------ est= uround ests= uround do 300 i=1,n s=e1*ak1(i)+e2*ak2(i)+e3*ak3(i)+e4*ak4(i) sk= dmax1(1.d-30,dabs(y(i)),dabs(ynew(i))) if (i.le.ny) est=dmax1(est,dabs(s)/sk) if (i.gt.ny) ests=dmax1(ests,dabs(s)/sk) 300 continue c---------------------new step size---------------------------------- tps= tol/est if (np.gt.0) tps=dmin1(tps,tols/ests) if (.not.end) hnew= 1dmin1(hmax,h*dmin1(facmx,dmax1(facmn,.9d0*tps**.25d0))) c-----------------------error test---------------------------------- if (end) goto 51 if (est.gt.tol) goto 43 if ( (np.gt.0).and.(ests.gt.tols) ) goto 43 c-------------------------accept step----------------------------- if (t.gt.tend) goto 48 facmn= .5d0 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,ests 100 format(5x,i5,' t= ',1pd10.3,' h= ',1pd10.3, 1 ' est= ',1pd10.3,' ests= ',1pd10.3) if(t.eq.tend) return goto 11 c-----------------------test for back step------------------------- 43 if (ipas.gt.0) irep=irep+1 if (hnew.gt.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-------------------------- 47 h= hnew facmx= 1.d0 goto 11 c--------------------manage the final step------------------------- 48 h= tend-told end= .true. goto 44 c--------------------succesful 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 * LINS sets up linear equations * c * * c ***************************************************************** c c np1= np + 1 c c common icue(np1),ipx(nec),ipy(nec),e2(nec) c subroutine lins(np,ny,h,y,fy,e,f2yp,inp) implicit real*8(a-h,o-z) implicit integer*2 (i-n) dimension y(ny,1),fy(ny,1),e(ny,1),f2yp(ny,1),inp(1) common/eec/ icue(21),ipx(2000),ipy(2000),e2(2000) c nec= 2000 ndy= ny if(np.eq.0) call jfy(ndy,y,fy) s= .395*h do 2 i=1,ny do 1 j=1,ny 1 e(i,j)=-s*fy(i,j) 2 e(i,i)= 1.d0+e(i,i) if (np.eq.0) return icue(1)= 0 do 3 k=1,np k1= k+1 call jf2yp(ndy,y,f2yp,inp(k)) do 4 i=1,ny call jf2yy(ndy,y,fy,i) do 4 j=1,ny d= f2yp(i,j) do 5 l=1,ny 5 d= d + fy(j,l) * y(l,k1) 4 f2yp(i,j)= s * d ic= icue(k) do 6 i=1,ny do 6 j=1,ny if(dabs(f2yp(i,j)).lt.1.d-30) goto 6 ic= ic + 1 if(ic.gt.nec) goto 7 ipx(ic)= i ipy(ic)= j e2 (ic)= f2yp(i,j) 6 continue icue(k+1)= ic 3 continue c write(6,*) ' ic = ',ic return 7 write(6,100) nec 100 format(//2X,'ERROR in subroutine LINS '// 1 2x,' EEC common must be redimensioned '// 2 2x,' in LINS and SOLS subroutines '// 3 2x,' Be new dimension greater then',i4,' !'// 4 2x,' --- PROGRAM TERMINATED --- '///) stop end c c *** KINAL ******************************************************* c * * c * FCNS computes the extended function * c * * c ***************************************************************** c subroutine fcns(np,ny,y,dy,fy,inp) implicit real*8(a-h,o-z) implicit integer*2(i-n) dimension y(ny,1),dy(ny,1),fy(ny,1),inp(1) c call fun(y,dy) if (np.eq.0) return ndy= ny call jfy(ndy,y,fy) call jfp(np,ndy,y,dy(1,2),inp) np1= np+1 do 1 k=2,np1 1 call au(ny,dy(1,k),fy,y(1,k)) return end c c *** KINAL ****************************************************** c * * c * SOLS solves extended system of linear equations * c * * c **************************************************************** c subroutine sols(np,ny,e,ee,ip,ak) implicit real*8(a-h,o-z) implicit integer*2 (i-n) dimension e(ny,1),ee(ny,1),ip(1),ak(ny,1) common/eec/ icue(21),ipx(2000),ipy(2000),e2(2000) c call sol(ny,e,ip,ak) if (np.eq.0) return do 1 k=1,np k1= k+1 do 2 i=1,ny do 2 j=1,ny 2 ee(i,j)= 0. do 3 i=(icue(k)+1),icue(k+1) 3 ee(ipx(i),ipy(i))=e2(i) call au(ny,ak(1,k1),ee,ak) 1 call sol(ny,e,ip,ak(1,k1)) return end c c *** KINAL ****************************************************** c * * c * AU is an auxiliary subroutine * c * * c **************************************************************** c subroutine au(n,a,b,c) implicit real*8(a-h,o-z) implicit integer*2 (i-n) dimension a(1),b(n,1),c(1) do 2 i=1,n d= a(i) do 1 j=1,n 1 d= d+b(i,j)*c(j) 2 a(i)= d 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 * 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 * 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 * JFP computes the Jacobian with respect to parameters * c * * c **************************************************************** c c y - IN: actual concentrations c fp - OUT: Jacobian matrix c fp(i,k) = d f(i) / d p(k) c ndy - IN: leading dimension of matrix fp c inp - IN: serial number of parameters , c related to which sensitivities are to be computed 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 jfp (np,ndy,y,fp,inp) 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),fp(ndy,1),r(90),inp(1) c do 1 k=1,nr 1 r(k)= 1. 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 do 3 k=1,np 3 fp(i,k)= 0. do 4 is=1,nisa do 4 k=1,np 4 if(inp(k).eq.isa(1,is)) 1 fp(isa(2,is),k)= - dble(isa(3,is)) * r(isa(1,is)) do 5 is=1,nisb do 5 k=1,np 5 if(inp(k).eq.isb(1,is)) 1 fp(isb(2,is),k)= sb(is) * r(isb(1,is)) return end c c *** KINAL ****************************************************** c * * c * JF2YP computes the mixed second derivatives * c * * c **************************************************************** c c y - IN: actual concentrations c f2yp - OUT: k1-th layer of the 3D array of second derivatives c f2yp(i,j) = d2 f(i) / d y(j) d p(k1) c k1 - IN: index of the layer (see above) c ndy - IN: leading dimension of matrix f2yp c q - working area c c ---------------------------------------------------------------- c c dimension q(ny) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine jf2yp(ndy,y,f2yp,k1) 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 f2yp(ndy,1),y(1),q(50) c do 1 i=1,ny do 1 j=1,ny 1 f2yp(i,j)=0. if(k1.gt.nr) return do 2 i=1,ny 2 q(i) = 0. do 3 is=1,nisa 3 if(isa(1,is).eq.k1) q(isa(2,is))= 1. do 4 j=1,ny do 4 is=1,nisa if ( isa(1,is) .ne.k1) goto 4 if ( isa(2,is) .eq. j) then q(j)= q(j) * real(isa(3,is)) * y(isa(2,is)) ** (isa(3,is)-1) else q(j)= q(j) * y(isa(2,is)) ** isa(3,is) endif 4 continue do 6 j=1,ny do 5 is=1,nisa if ( isa(1,is) .ne.k1) goto 5 f2yp(isa(2,is),j)= f2yp(isa(2,is),j) - real(isa(3,is)) * q(j) 5 continue do 6 is=1,nisb if ( isb(1,is) .ne.k1) goto 6 f2yp(isb(2,is),j)= f2yp(isb(2,is),j) + sb(is) * q(j) 6 continue return end c c *** KINAL ****************************************************** c * * c * JF2YY computes the second derivatives * c * with respect to variables * c * * c **************************************************************** c c y - IN: actual concentrations c f2yy - OUT: i1-th layer of the 3D array of second derivatives c f2yy(i,j) = d2 f(i1) / d y(i) d y(j) c i1 - IN: index of layer (see above) c ndy - IN: leading dimension of matrix f2yy c dr,ir,iq - working areas c c ---------------------------------------------------------------- c c dimension dr(nr),ir(nr),iq(nr) c common isa(3,nca),isb(2,ncb),sb(ncb),p(nr) c subroutine jf2yy(ndy,y,f2yy,i1) 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 f2yy(ndy,1),y(1),dr(90),ir(90),iq(90) c do 1 k=1,nr 1 iq(k)= 0 do 2 is=1,nisa 2 if(isa(2,is).eq.i1) iq(isa(1,is))= 1 do 3 is=1,nisb 3 if(isb(2,is).eq.i1) iq(isb(1,is))= 1 c do 14 i=1,ny do 14 j=1,i do 4 k=1,nr ir(k)= 0 4 dr(k)= 0. f2yy(i,j)= 0. if(i.eq.j) goto 7 c -------------- off diagonal elements --------------------- do 5 is=1,nisa if (iq(isa(1,is)).ne.1) goto 5 ip= (isa(2,is)-j) * (isa(2,is)-i) in1= isa(1,is) if (ip.eq.0) ir(in1)= ir(in1) + 1 if (ir(in1).ge.2) dr(in1)= p(in1) 5 continue do 6 is=1,nisa if (iq(isa(1,is)).ne.1) goto 6 if (ir(isa(1,is)).lt.2) goto 6 ip= (isa(2,is)-j) * (isa(2,is)-i) if (ip.eq.0) then dr(isa(1,is)) = dr(isa(1,is))*real(isa(3,is)) 1 *y(isa(2,is))**(isa(3,is)-1) else dr(isa(1,is)) = dr(isa(1,is))*y(isa(2,is))**isa(3,is) endif 6 continue goto 10 c -------------- diagonal elements ------------------------- 7 do 8 is=1,nisa if (iq(isa(1,is)).ne.1) goto 8 if (isa(2,is).ne.j) goto 8 in1= isa(1,is) if(isa(3,is).ge.2) then ir(in1)= 2 dr(in1) = p(in1) endif 8 continue do 9 is=1,nisa if (iq(isa(1,is)).ne.1) goto 9 if (ir(isa(1,is)).lt.2) goto 9 if (isa(2,is).eq.j) then dm= dble(isa(3,is)) dr(isa(1,is)) = dr(isa(1,is))*dm*(dm-1.)* 1 y(isa(2,is))**(isa(3,is)-2) else dr(isa(1,is)) = dr(isa(1,is))*y(isa(2,is))**isa(3,is) endif 9 continue c -------------- summation --------------------------------- 10 continue do 11 is=1,nisa if(isa(2,is).ne.i1) goto 11 if(dr(isa(1,is)).eq.0.) goto 11 f2yy(i,j)= f2yy(i,j) - real(isa(3,is)) * dr(isa(1,is)) 11 continue do 12 is=1,nisb if(isb(2,is).ne.i1) goto 12 if(dr(isb(1,is)).eq.0.) goto 12 f2yy(i,j)= f2yy(i,j) + sb(is) * dr(isb(1,is)) 12 continue 14 f2yy(j,i)= f2yy(i,j) 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