c c ***************************************************************** c * * c * MECHMOD version 3.42 * c * * c * Utility program to the CHEMKIN-III package * c * for the transformation of mechanism files * c * * c ***************************************************************** c c The latest version of this program is always available from c the World Wide Web: c http://www.chem.leeds.ac.uk/Combustion/Combustion.html c or c http://garfield.chem.elte.hu/Combustion/Combustion.html c c---------------------------------------------------------------------- c c This program has to be linked to the CHEMKIN_LIBRARY c of the CHEMKIN-III (version 3.6) package. c Please use makefile mechmod4CKIII.mak c c---------------------------------------------------------------------- c c Version 3.42 (21st February, 2003) c c Changes from version 3.41 to version 3.42 (21st February, 2003) c 1 Multiple lines of 3rd body coefficients are handled c c Changes from version 3.40 to version 3.41 (16th May, 2002) c 1 The chemical interpreter of CHEMKIN-III checks duplicated c third body reactions (unlike that of CHEMKIN-II). c MECHMOD for CHEMKIN-III was modified accordingly. c c Changes from version 1.4 to version 3.40 (5th April, 2002) c 1 MECHMOD version 1.4 that uses the CHEMKIN-II subroutine cklib.f c was modified to work with the CHEMKIN-III package. c In this version only CHEMKIN-II compatible mechanisms are transformed. c New options of the mechanism files introduced in CHEMKIN-III c are not interpreted yet. c c Changes from version 1.3 to version 1.4 (31st March, 2001) c 1 New keywords THINFO_kJ, THINFO_kcal c for printing thermodynamic info summary in the mechanism file c 2 New keywords to modify the 298.15K thermodynamic data c stored in the NASA polynomials: c newH_kJ, newH_kcal, newS_J/K, newS_cal/K c 3 Command KILL is accepted multiple times c 4 Species given in KILL are eliminited from c the list of species and from the list of 3rd body species c 5 Initial value of lWL is set to .false. c 6 Handling of Troe parameters in normal form c 7 NoREA was added to the argument list of subroutine rwrite c 8 Deleting 'lout' from the argument list of backA and backLT c c Changes from version 1.2 to version 1.3 (20th February, 1998) c 1 LOUT was added to the argument list of subroutine rwrite c 2 An interference between keywords KILL and THERMO was eliminated c 3 Parameter n is expressed with more significant digits. c 4 More precise gas constant is applied, c if the CKINTERP used is v.3.8 (6/13/94) or later. c 5 COMMON CKSTRT in subroutine strt1 was extended with a dummy c array to avoid warning messages. c 6 3rd body efficiency string is printed in a shorter form. c c Changes from version 1.1 to version 1.2 (12th July, 1995) c 1 The structure of COMMON CKSTRT in the CHEMKIN package was c changed on 1/26/94 (from CKINTERP v.3.3 and CKLIB v.4.5). c To keep the compatibility, a version sensitive handling c of CKSTRT has been introduced to MECHMOD. c 2 More space is allowed for printing the temperature exponents. c c Changes from version 1.0 to version 1.1 (8th June, 1995) c 1 The program reads the version of the CKINTERP program that c was used for the creation of the cklink file and takes into c account the version number. c This allows utilization of CHEMKIN versions prior to 15/2/91 c (prior to CKINTERP ver. 2.6 and CHEMKIN ver. 2.9) c c Version 1.0 was launched on 1st May, 1995 c c------------------------------------------------------------------------ c c The program creates a mechanism text file from a CHEMKIN-III c chem.asc file. If the options below are not used, c the new mechanism file differs from the original as follows: c - the comments are missing c - the format may be nicer c - the units become mole-cm-sec-K for the prexponential factors and c degrees Kelvin for the activation energies. c c Options: c c - All reversible reactions are transformed to pairs of c irreversible reactions c c - Selected species are eliminated from the mechanism c c - The preexponential factors and the activation energies c are converted to different units c c - The thermodynamic data, stored in the CHEMKIN chem.bin file, c are printed at the beginning of the mechanism c in the form of NASA polynomials c c - Modification of room temperature (298.15K) thermodynamic data c c - Thermodynamic data summary is printed in the mechanism c as comment lines c c c Subroutine CKSYMRL in this program is a variant of subroutine CKSYMR c of the CHEMKIN package (ver. 4.9). Subroutines UPCAS (originally c UPCASE) and CKDUP were loaned from program CKINTERP version 3.6. c c Version 1.0 of this program was developed using CHEMKIN-II version 4.9 c c************************************************************************ c c Please feel free to distribute this program c under the following conditions: c c - Distribution of modified versions is not allowed c - All bugs found have to be reported to c the author: c c Tamas Turanyi c c Department of Physical Chemistry c Eotvos University (ELTE) c H-1518 Budapest, P.O.Box 32, Hungary c Phone : (36-1) 209-0555 x1109 c Fax : (36-1) 209-0602 c E-mail : turanyi@garfield.chem.elte.hu c turanyi@chemres.hu c tamast@chem.leeds.ac.uk c WWW : http://garfield.chem.elte.hu/ c c----------------------------------------------------------------------- c c Matching to the CHEMKIN-III package was written by c c Istvan Gy. Zsely c c Department of Physical Chemistry c Eotvos University (ELTE) c E-mail : zsigy@vuk.chem.elte.hu c WWW : http://garfield.chem.elte.hu/ c c------------------------------------------------------------------- c c The author is very grateful for the helpful suggestions to c c Jean-Claude Boettner (Orleans) c Melita Jazbec (Sydney) c Zsely, Istvan Gyula (Budapest) c Pankaj Bajaj (Zurich) c Kevin Hughes (Leeds) c c=================================================================== c c c C O N T R O L c c c The program accepts the keywords below, which have to be c in a text file (called 'control file'). c c c c> IRREV transforms all the reversible reactions to pairs of c irreversible reactions. c c This keyword may be followed by an integer number 'n', c which controls the way that the reverse Arrhenius parameters c are calculated. c c n=1 Reverse Arrhenius parameters A, n, and E are fitted c at temperatures Tlow, Tmid, and Thigh. These temperatures c are defined in the first executable lines of the program. c This is the default. c c n=2 Forward Arrhenius parameter n is preserved as c reverse parameter n. Reverse Arrhenius parameters A and E c are fitted at temperatures Tlow and Thigh. c c n=3 Forward Arrhenius parameter A is preserved as c reverse parameter A. Reverse Arrhenius parameters n and E c are fitted at temperatures Tlow and Thigh. c c> THERMO All thermodynamic data are printed at the beginning of the c mechanism. The default is that the thermodynamic data are c not given at the beginning of the mechanism. c c> KILL Eliminates those reactions from the mechanism c that include (as reactants or products) the named species. c Followed by the list of species to be eliminated. c c> FORMAT1 The reactions are printed as follows: c CH4 + OH => CH3 + H2O 1.570E07 1.8300 11.64 c This is the default c c> FORMAT2 The reactions are printed as follows: c CH4 + OH => CH3 + H2O 1.570E07 1.830 11.64 c c> MOLES The printed preexponential units are moles-cm-sec-K (default) c c> MOLECULES The printed preexponential units are molecules-cm-sec-K c c> KELVINS The printed activation energies are Kelvin (default) c c> CAL/MOLE The printed activation energies are cal/mole c c> KCAL/MOLE The printed activation energies are kcal/mole c c> JOULES/MOLE The printed activation energies are Joules/mole c c> KJOULES/MOLE The printed activation energies are kJoules/mole. c c> THINFO_kcal Heat-of-formation and entropy of species at 298.15K c are printed in the mechanism file c in kcal and cal/K units, respectively. c c> THINFO_kJ Heat-of-formation and entropy of species at 298.15K c are printed in the mechanism file c in kJ and J/K units, respectively. c c c The following commands modify the thermodynamic data stored in NASA polynomials. c The new heat-of-formation and entropy data refer to 298.15K. c c> newH_kJ c c> newH_kcal c c> newS_J/K c c> newS_cal/K c c> END The control file may be closed by an END statement. c The subsequent commands are ignored. c c c Further rules: c c - the order of keywords is arbitrary (except for END) c - one keyword is allowed on each line c - characters after an '!' sign are not interpreted c - lower case and capital letters are not distinguished c c An example for a valid control file: c c IRREV 3 c kill H2O CO CO2 c !FORMAT2 c MOLECULES c KJOULES/MOLE ! unit for E is kJ/mole c THINFO_KJ c newH_kJ OH 41.0 c END c c------------------------------------------------------------------------ c c linc = Unit number for CHEMKIN linking file input c lctr = Unit number for the control keywords input c lkb = Unit number for keyboard input c lout = Unit number for screen messages c lmech = Unit number for the created mechanism c lnul = Unit number for discarded error messages c c leni = Length of integer work array c lenr = Length of real work array c lenc = Length of character work array c lensym = Length of a character string in character work array c c iw = Integer work array c rw = Real work array c cw = Character work array c c logical variables: c c lrev .true. reversible reactions are permitted in the output mech. c .false. all reversible reactions are converted to irrev. c c lirrev .true. the reaction is originally irreversible c .false. the reaction is originally reversible c c ltherm .true. the thermodynamic data are printed at the c beginning of the mechanism c c L16 .true. at least one of the species' names are c more than 7 characters long c c kerr .true. error in a CHEMKIN subroutine c c lArev .true. reverse Arrhenius parameters were given to c a reversible reaction c c lLT .true. Landau-Teller parameters were given to the reaction c c lLrev .true. reverse Landau-Teller parameters were given to c a reversible reaction c c l3rd .true. enhanced 3rd body info for the reaction c c ldup .true. the reaction is duplicate c c lkill .true. the reaction contains an unwanted species c c lWL .true. reaction with radiation wavelength c c------------------------------------------------------------------------- c implicit double precision (a-h, o-z), integer (i-n) parameter ( leni =18000, lenr =15000, lenc = 300, 1 lensym = 16) parameter (mxkill=100, mxdup=100, mxmod=100) dimension iw(leni), rw(lenr) character cw(lenc)*(lensym), cspec(lenc)*(lensym) dimension ikill(mxkill),idup(mxdup), iHmod(mxmod), iSmod(mxmod) dimension rHmod(mxmod), rSmod(mxmod), aup(7),alp(7),ke(5) character fdflt*10, che(5)*2, chph*1, vers*16 character*80 fname, line, istrl, istrr, lstr character gline(10)*80 logical kerr,L16,lrev,ltherm,lirrev, 1 lArev,lLT,lLrev,l3rd,ldup,lWL,lkill data linc/3/, lctr/4/, lkb/5/, lout/6/, lmech/7/, lnul/8/ c include 'ckstrt.h' c c Reverse parameters are fitted at temperatures c Tlow, Tmid, and Thigh [K] c Tlow= 1000. Tmid= 1750. Thigh= 2500. c c Opening I/O channels c write(lout,200) 200 format(// 1 5x,'MECHMOD'/ 2 5x,'Program for the modification of ', 3 'CHEMKIN-III mechanism files'/ 4 5x,'Version 3.42 (21st February, 2003)'// 5 5x,'Hit ENTER for defaults:'//) c fdflt='chem.asc' write(lout,201) fdflt 201 format(1x,'Name of CHEMKIN-III link file [',a10,'] :') read (lkb ,100) fname if (fname.eq.' ') fname= fdflt open (unit = linc, status= 'old', form= 'formatted', 1 file= fname, iostat= mes, err= 410) c fdflt= 'Mechanism' write(lout,202) fdflt 202 format(1x,'Name of the file of the new mechanism [',a10,'] :') read (lkb ,100) fname if (fname.eq.' ') fname= fdflt open (unit = lmech, status= 'unknown', form= 'formatted', 1 file= fname, iostat= mes, err= 420) c fdflt= 'mcontrol' write(lout,203) fdflt 203 format(1x,'Name of the file of keywords [',a10,'] :') read (lkb ,100) fname if (fname.eq.' ') fname= fdflt open (unit = lctr, status= 'old', form= 'formatted', 1 file= fname, iostat= mes, err= 430) c c UNIX c c open (unit = lnul, status='old', file='/dev/null') c c DOS / Windows c open (unit = lnul, status='old', file='NUL') c c--------------------------------------------------------------------- c c determine mechanism size and initialization of CHEMKIN c call cklen (linc, lout, leniwk, lenrwk, lencwk, iflag) c if (leniwk.le.leni.and.lenrwk.le.lenr.and.lencwk.le.lenc) then call ckinit (leniwk,lenrwk,lencwk,linc,lout,iw,rw,cw,iflag) else write (lout,204) leni, leniwk, lenr, lenrwk, lenc, lencwk write (lout,205) stop endif 204 format (/, ' Working Space Requirements', 1 /, ' Provided Required ', 2 /, ' Integer ', 2I15, 3 /, ' Real ', 2I15, 4 /, ' Character', 2I15, /) 205 format (' Stop, not enough work space provided.') c zsigy NPAR4=NPAR+1 c--------------------------------------------------------------------- c c max. length of species names c L16= .false. do 1 k=1,NKK ilen= ilasch(cw(IcKK+m-1)) if(ilen.ge.8) L16= .true. 1 continue c c-------------------------------------------------------------------- c call control(lctr,lout,lnul,NKK,cw(IcKK),L16,iformat,ltherm,lrev, 1 metrev,nkill,mxkill,ikill,iunitA,iunitE, 2 ithinfo, mxmod, nHmod, iHmod, rHmod, nSmod, iSmod, rSmod, 3 Tlow,Tmid,Thigh) c c---------------------------------------------------------------------- c c location of DUPLICATE reactions c ndup= 0 do 2 i=1,NII CALL CKDUP (i, MXSP, iw(IcNS), iw(IcNR), iw(IcNU), iw(IcNK), 1 iw(IcFL), iw(IcKF), iw, rw, isame) c if (isame.ne.0) then ndup= ndup+2 idup(ndup)= isame idup(ndup-1)= i c c A+B(+M) and A+B(+C) are not duplicates c i1= -1 i2= -1 do 25 n=1,NFAL if(iw(IcFL+n-1).eq.i) i1= iw(IcKF+n-1) if(iw(IcFL+n-1).eq.isame) i2= iw(IcKF+n-1) 25 continue if (i1.ne.i2) ndup= ndup-2 endif 2 continue c c---------------------------------------------------------------------- c c writing the head of the mechanism file c line= 'ELEMENTS' do 26 m=1,NMM ilen= ilasch(cw(IcMM+m-1)) ilenl= ilasch(line) line= line (:ilenl+2) // cw(IcMM+m-1)(:ilen+1) 26 continue ilenl=ilasch(line) line= line(:ilenl+2) // 'END' write(lmech,300) line 300 format(a80) c write(lmech,301) 301 format('SPECIES') m1= 0 do 36 m=1,NKK do 37 k=1,nkill 37 if (ikill(k).eq.m) goto 36 m1=m1+1 cspec(m1)= cw(IcKK+m-1) 36 continue c if (L16) then write(lmech,302)(cspec(m),m=1,m1) 302 format(5a16) else write(lmech,303)(cspec(m),m=1,m1) 303 format(10a8) endif write(lmech,304) 304 format('END') c c-------------------------------------------------------------------- c c writing thermodynamic data summary / modification of data c if (ithinfo.eq.0) goto 59 write(lmech,550) if (ithinfo.eq.1) write(lmech,551) if (ithinfo.eq.2) write(lmech,552) 550 format('!',/ 1 '! Heat-of-formation, entropy, and ', 2 'heat capacity of species at 298.15K'/'!'/ 3 '!'/'! species dH S ', 4 'Cp Tlow Tmid Thigh') 551 format( '! kJ/mole J/K_mole ', 1 'kJ/mole K K K'/'!') 552 format( '! kcal/mole cal/K_mole ', 1 'kcal/mole K K K'/'!') c do 50 k=1,NKK c c temperature range c TLO= rw(NcTT+(k-1)*MXTP ) TMI= rw(NcTT+(k-1)*MXTP+1) THI= rw(NcTT+(k-1)*MXTP+2) c do 51 i=1,7 alp(i)= rw(NcAA+ (k-1)*NCP2T+i-1) aup(i)= rw(NcAA+NCP2+(k-1)*NCP2T+i-1) 51 continue c T298=298.15 R=8.314 c c lower temperature range polynomial set c sumH = 0. sumS = 0. sumC = 0. c sumH = sumH + alp(6) do 52 i=1,5 sumH = sumH + alp(i)*(T298**i)/i sumC = sumC + alp(i)*(T298**(i-1)) 52 continue C sumS = alp(1)*dlog(T298) + alp(7) do 53 i=2,5 sumS = sumS + alp(i)*T298**(i-1)/(i-1) 53 continue H1 = sumH*R/1000. S1 = sumS*R C1 = sumC*R c c upper temperature range polynomial set c sumH = 0. sumS = 0. sumC = 0. c sumH = sumH + aup(6) do 54 i=1,5 sumH = sumH + aup(i)*(T298**i)/i sumC = sumC + aup(i)*(T298**(i-1)) 54 continue C sumS = aup(1)*dlog(T298) + aup(7) do 55 i=2,5 sumS = sumS + aup(i)*T298**(i-1)/(i-1) 55 continue c H2 = sumH*R/1000. S2 = sumS*R C2 = sumC*R c c modification of thermodynamic data c if (nHmod.eq.0) goto 60 do 56 i=1,nHmod if(iHmod(i).eq.k) then alp(6)=alp(6)+(rHmod(i)-H1)/R*1000 aup(6)=aup(6)+(rHmod(i)-H1)/R*1000 H1o=H1 H1=rHmod(i) H2=rHmod(i)+H2-H1o endif 56 continue 60 continue if (nSmod.eq.0) goto 61 do 57 i=1,nSmod if(iSmod(i).eq.k) then alp(7)=alp(7)+(rSmod(i)-S1)/R aup(7)=aup(7)+(rSmod(i)-S1)/R S1o=S1 S1=rSmod(i) S2=rSmod(i)+S2-S1o endif 57 continue 61 continue c do 62 i=1,7 rw(NcAA+ (k-1)*NCP2T+i-1) = alp(i) rw(NcAA+NCP2+(k-1)*NCP2T+i-1) = aup(i) 62 continue c c printing c if (ithinfo.eq.1) then write(lmech,553) cw(IcKK+k-1),H1,S1,C1,TLO,TMI,THI c write(lmech,553) cw(IcKK+k-1),H2,S2,C2,TLO,TMI,THI else c2J= 4.184 write(lmech,553) cw(IcKK+k-1), H1/c2J, S1/c2J, C1/c2J, 1 TLO,TMI,THI c write(lmech,553) cw(IcKK+k-1), H2/c2J, S2/c2J, C2/c2J, c 1 TLO,TMI,THI endif 553 format('! ',a16,3f11.3,3x,3f7.0) c 50 continue write(lmech,554) 554 format('! '/ 1 '! Thermodynamic summary was printed by program', 2 ' MECHMOD v3.42'/'!') 59 continue c c----------------------------------------------------------------------c c the thermodynamic data c if (ltherm) then write(lmech,305) 305 format('THERMO ALL') TLO= 300.0 TMI= 1000.0 THI= 5000.0 write(lmech,306) TLO,TMI,THI 306 format(3f10.3) c c identification of elemental composition do 40 k=1,NKK do 41 i=1,5 che(i)= ' ' 41 ke(i)= 0 ie= 0 do 42 m=1,NMM kes= iw(IcNC+(k-1)*NMM+M-1) if (kes.ne.0) then ie=ie+1 che(ie)= cw(IcMM+m-1)(:2) ke (ie)= kes endif 42 continue c if (iw(IcCH+k-1).ne.0) then ie=ie+1 che(ie)= 'E ' ke (ie)= iw(IcCH+k-1) endif c chph= 'G' if (iw(IcPH+k-1).eq.-1) chph= 'S' if (iw(IcPH+k-1).eq. 1) chph= 'L' c TLO= rw(NcTT+(k-1)*MXTP ) TMI= rw(NcTT+(k-1)*MXTP+1) THI= rw(NcTT+(k-1)*MXTP+2) c do 43 i=1,7 alp(i)= rw(NcAA+ (k-1)*NCP2T+i-1) aup(i)= rw(NcAA+NCP2+(k-1)*NCP2T+i-1) 43 continue c write(lmech,307) cw(IcKK+k-1),(che(i),ke(i),i=1,4),chph, 1 TLO,THI,TMI,che(5),ke(5),1 307 format(a16,8x,4(a2,i3),a1,2f10.2,f8.2,a2,i3,1x,i1) write(lmech,308) (aup(i),i=1,5),2 write(lmech,308) (aup(i),i=6,7),(alp(i),i=1,3),3 write(lmech,309) (alp(i),i=4,7),4 308 format(5(1pe15.8),4x, i1) 309 format(4(1pe15.8),19x,i1) 40 continue write(lmech,304) endif c c the reactions c NoREA= 0 line= 'REACTIONS ' if (iunitA.eq.1) line=line(:12)//'MOLES ' if (iunitA.eq.2) line=line(:12)//'MOLECULES ' ilenl= ilasch(line) if (iunitE.eq.1) line=line(:ilenl)//' KELVINS' if (iunitE.eq.2) line=line(:ilenl)//' CAL/MOLE' if (iunitE.eq.3) line=line(:ilenl)//' KCAL/MOLE' if (iunitE.eq.4) line=line(:ilenl)//' JOULES/MOLE' if (iunitE.eq.5) line=line(:ilenl)//' KJOULES/MOLE' c write(lmech,300) line do 3 i=1,NII c c--------------------------------------------------------------- c Setting some reaction status flags c c lirrev= .true. the reaction is irreversible c in the original mechanism c lirrev= .false. if (iw(IcNS+i-1).lt.0) lirrev= .true. c c l3rd= .true. there is a 3rd body in the reaction c l3rd= .false. do 8 n=1,NTHB if (i.eq.iw(IcTB+n-1)) then l3rd= .true. n3rd= n endif 8 continue c c ifltyp.ne.0 this is a fall-off reaction c ifltyp= 0 do 19 n=1,NFAL if (i.eq.iw(IcFL+n-1)) then ifltyp= iw(IcFO+n-1) nfl= n endif 19 continue c c ldup= .true. this is a duplicate reaction c ldup= .false. do 28 n=1,ndup 28 if(i.eq.idup(n)) ldup= .true. c c lkill= .true. this reaction will be eliminated c lkill= .false. do 29 k=1,nkill do 29 n=1,MXSP 29 if(iw(IcNK+(i-1)*MXSP+n-1).eq.ikill(k)) lkill= .true. c c--------------------------------------------------------------- c c the basic reaction line c A= rw(NcCO+(i-1)*NPAR4+0) Tn= rw(NcCO+(i-1)*NPAR4+1) E= rw(NcCO+(i-1)*NPAR4+2) ico= 1 call convert(i,MXSP,iw(IcNU),ico,ifltyp,l3rd, 1 iunitA,iunitE,A,E,Ac,Ec,rvers) c JR= 1 call CKSYMRL (I, LOUT, IW, RW, CW, JR, LTL, ISTRL, RVERS, KERR) if(kerr) stop JR= 2 call CKSYMRL (I, LOUT, IW, RW, CW, JR, LTR, ISTRR, RVERS, KERR) if(kerr) stop c c reaction output c if(lkill) then line= '! Deleted reaction: '//ISTRL(:LTL)//' ='//ISTRR(:LTR) write(lmech,325) 325 format('!') write(lmech,300) line goto 3 else NoREA= NoREA+1 write(lmech,319) NoREA 319 format('!'/'!',I6) call rwrite(lout,lmech,iformat,lrev,lirrev,istrl,istrr, 1 Ac,Tn,Ec,iferr,NoREA) if (iferr.eq.1) then call rwrite(lout,lmech,iformat,lrev,lirrev,istrl,istrr, 1 Ac,Tn,Ec,iferr,NoREA) iformat=2 endif endif c c duplicate reactions c if(ldup) write(lmech,318) 318 format('DUPLICATE') c c reverse Arrhenius parameters c lArev= .false. do 4 n=1,NREV if (i.eq.iw(IcRV+n-1)) then c reverse Arrhenius parameters lArev= .true. Ar= rw(NcRV+(n-1)*NPAR4+0) Tnr= rw(NcRV+(n-1)*NPAR4+1) Er= rw(NcRV+(n-1)*NPAR4+2) endif 4 continue if(lrev .and. lArev) then ico= 2 call convert(i,MXSP,iw(IcNU),ico,ifltyp,l3rd, 1 iunitA,iunitE,Ar,Er,Arc,Erc,rvers) write(lmech,310) Arc,Tnr,Erc 310 format(' REV / ',1pe10.3,0pf6.2,f10.2,' /') endif c c Landau-Teller reactions c lLT= .false. do 5 n=1,NLAN if (i.eq.iw(IcLT+n-1)) then c LT parameters lLT= .true. TL1= rw(NcLT+(n-1)*NLAR+0) TL2= rw(NcLT+(n-1)*NLAR+1) write(lmech,311) TL1,TL2 311 format(' LT / ',2(1pe10.3),' /') endif 5 continue c c Reverse LT data c lLrev= .false. do 6 n=1,NRLT if (i.eq.iw(IcRL+n-1)) then c reverse LT parameters lLrev= .true. TL1r= rw(NcRL+(n-1)*NLAR+0) TL2r= rw(NcRL+(n-1)*NLAR+1) endif 6 continue if(lrev.and.lLrev) write(lmech,312) TL1r,TL2r 312 format(' RLT / ',2(1pe10.3),' /') c c wavelength information c lWL= .false. do 7 n=1,NWL if (i.eq.iw(IcWL+n-1)) then c wavelength information lWL= .true. WL= rw(NcWL+n-1) write(lmech,313) -WL 313 format(' HV / ',1pg10.3,' /') endif 7 continue c c enhanced 3rd bodies c do 1000 l1=1,10 gline(l1)= ' ' 1000 continue if (.not.l3rd) goto 91 if (iw(IcKN+n3rd-1).eq.0) goto 91 c if l1=1 then this is the first 3rd body in the line l1=0 nline=1 do 9 l=1,iw(IcKN+n3rd-1) l1=l1+1 i3rdsp= iw(IcKT+(n3rd-1)*MXTB+l-1) c c eliminating unwanted 3rd body species c do 35 k=1,nkill if (ikill(k).eq.i3rdsp) goto 9 35 continue c ilenc= ilasch(cw(IcKK+i3rdsp-1)) r3rd= rw(NcKT+(n3rd-1)*MXTB+l-1) call r2ch(r3rd,lstr,ilenr,kerr) if (kerr) stop ilenl= ilasch(gline(nline))+1 if ((ilenl+ilenc+7).gt.80) then write(lmech,300) gline(nline) nline=nline+1 l1=1 endif if (l1.eq.1) ilenl=1 gline(nline)= gline(nline)(:ilenl) 1 // cw(IcKK+i3rdsp-1)(:ilenc) // '/'// lstr(:ilenr) // '/' 9 continue write(lmech,300) gline(nline) 91 continue c c fall-off reactions c if (ifltyp.eq.0) goto 10 c c parameters of the fall-off reactions: c Al= rw(NcFL+(nfl-1)*NFAR+0) Tnl= rw(NcFL+(nfl-1)*NFAR+1) El= rw(NcFL+(nfl-1)*NFAR+2) p4= rw(NcFL+(nfl-1)*NFAR+3) p5= rw(NcFL+(nfl-1)*NFAR+4) p6= rw(NcFL+(nfl-1)*NFAR+5) p7= rw(NcFL+(nfl-1)*NFAR+6) p8= rw(NcFL+(nfl-1)*NFAR+7) c ico= 3 call convert(i,MXSP,iw(IcNU),ico,ifltyp,l3rd, 1 iunitA,iunitE,Al,El,Alc,Elc,rvers) write(lmech,314) Alc,Tnl,Elc 314 format(' LOW / ',1pe10.3,1x,0pf6.2,1x,f10.2,' /') c c branching according to the type of fall-off parameters c goto (14,11,12,13) ifltyp goto 14 c c SRI c 11 write(lmech,315) p4,p5,p6,p7,p8 315 format(' SRI / ',5(1pe12.4),' /') goto 14 c c TROE with 6 parameters c 12 write(lmech,316) p4,p5,p6 316 format(' TROE / ',3(1pe12.4),' /') goto 14 c c TROE with 7 parameters c 13 write(lmech,317) p4,p5,p6,p7 317 format(' TROE / ',4(1pe12.4),' /') 14 continue 10 continue c c Writing the (forward) rate info is complete c--------------------------------------------------------------------- c Optional writing of the backward reactions c if(lrev.or.lirrev) goto 3 NoREA= NoREA+1 if (.not.lArev) then call backA(metrev,i,iw,rw,rw(NcK1),Tlow,Tmid,Thigh, 1 A,Tn,E,Ar,Tnr,Er,rvers,NKK) ico= 2 call convert(i,MXSP,iw(IcNU),ico,ifltyp,l3rd, 1 iunitA,iunitE,Ar,Er,Arc,Erc,rvers) write(lmech,326) NoREA 326 format('!'/'!',I6,' R') else write(lmech,326) NoREA endif c c reaction output c call rwrite(lout,lmech,iformat,lrev,lirrev,istrr,istrl, 1 Arc,Tnr,Erc,iferr,NoREA) if (iferr.eq.1) then call rwrite(lout,lmech,iformat,lrev,lirrev,istrr,istrl, 1 Arc,Tnr,Erc,iferr,NoREA) iformat=2 endif c c duplicate reactions c if(ldup) write(lmech,318) c c Landau-Teller reactions c if (.not.lLT) goto 15 if (.not.lLrev) call backLT (i,iw,rw,rw(NcK1),Tlow,Thigh, 1 A,Tn,E,TL1,TL2,Ar,Tnr,Er,TL1r,TL2r,rvers,NKK) write(lmech,311) TL1r,TL2r 15 continue c c wavelength information c if (.not.lWL) goto 16 write(lmech,313) WL 16 continue c c enhanced 3rd bodies c if (.not.l3rd) goto 17 write(lmech,300) line 17 continue c c fall-off reactions c if (ifltyp.eq.0) goto 3 call backA(metrev,i,iw,rw,rw(NcK1),Tlow,Tmid,Thigh, 1 Al,Tnl,El,Alr,Tnlr,Elr,rvers,NKK) ico= 4 call convert(i,MXSP,iw(IcNU),ico,ifltyp,l3rd, 1 iunitA,iunitE,Alr,Elr,Alrc,Elrc,rvers) write(lmech,314) Alrc,Tnlr,Elrc c c the type of fall-off parameters c goto (23,20,21,22) ifltyp goto 23 c c SRI c 20 write(lmech,315) p4,p5,p6,p7,p8 goto 23 c c TROE with 6 parameters c 21 write(lmech,316) p4,p5,p6 goto 23 c c TROE with 7 parameters c 22 write(lmech,317) p4,p5,p6,p7 23 continue c c end of writing the reverse reactions c------------------------------------------------------------------- c 3 continue write(lmech,304) if(lrev) then write(lmech,320) else write(lmech,321) goto (31,32,33) metrev goto 39 31 write(lmech,322) Tlow,Tmid,Thigh goto 39 32 write(lmech,323) Tlow,Thigh goto 39 33 write(lmech,324) Tlow,Thigh goto 39 endif c 320 format('!-----------------------------------------------------!'/ 1 '! The mechanism was reprinted by program MECHMOD 3.42 !'/ 2 '!-----------------------------------------------------!') 321 format('!-----------------------------------------------------!'/ 1 '! Reverse Arrhenius parameters of reactions !'/ 2 '! denoted by "R" were calculated by program MECHMOD. !') 322 format('! Parameters A, n, and E were fitted at temperatures !'/ 1 '! ', f6.1,' K, ',f6.1,' K, and ',f6.1,' K',14x,' !'/ 3 '!-----------------------------------------------------!') 323 format('! Parameters A and E were fitted at temperatures !'/ 1 '! ', f6.1,' K and ',f6.1,' K',12x,' !'/ 3 '!-----------------------------------------------------!') 324 format('! Parameters n and E were fitted at temperatures !'/ 1 '! ', f6.1,' K and ',f6.1,' K',12x,' !'/ 3 '!-----------------------------------------------------!') 39 close (linc) close (lmech) close (lctr) stop 410 write (lout,299) mes,linc stop 420 write (lout,299) mes,lmech stop 430 write (lout,299) mes,lctr stop c 100 format(80a) 299 format(5x,'Error No',i4,' at opening I/O channel No',i3// 1 5x,' --- PROGRAM TERMINATED ---') c end C C----------------------------------------------------------------------C C SUBROUTINE CKSYMRL(I, LOUT, ICKWRK, RCKWRK, CCKWRK, J, LT, ISTR, 1 RVERS, KERR) C C START PROLOGUE C C SUBROUTINE CKSYMRL (I, ICKWRK, RCKWRK, CCKWRK, J, LT, ISTR, RVERS, KERR)* C Returns a character string which describes the left or right C hand side of the Ith reaction, C and the effective length of the character string. C C INPUT C I - Reaction index. C Data type - integer scalar C LOUT - Output unit for printed diagnostics. C Data type - integer scalar C ICKWRK - Array of integer workspace. C Data type - integer array C Dimension ICKWRK(*) at least LENIWK. C RCKWRK - Array of real work space. C Data type - real array C Dimension RCKWRK(*) at least LENRWK. C CCKWRK - Array of character work space. C Data type - CHARACTER*16 array C Dimension CCKWRK(*) at least LENCWK. C J = 1 the left hand side is given C = 2 the right hand side is given C C RVERS - The version of CHEMKIN C C OUTPUT C ISTR - Character string describing the Ith reaction. C Data type - CHARACTER*(*) C LT - Number of characters in the reaction description. C Data type - integer scalar C KERR - Error flag; character length error will result in C KERR=.TRUE. C Data type - logical C C END PROLOGUE C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) C*****END precision > double C*****precision > single C IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION ICKWRK(*),RCKWRK(*) CHARACTER CCKWRK(*)*(*), ISTR*(*), IDUM*80 LOGICAL KERR, IERR C COMMON /CMIN/ CKMIN C c include 'ckstrt.h' NPAR4=NPAR+1 C ISTR = ' ' ILEN = LEN(ISTR) KERR = .FALSE. C IRNU = 0 DO 10 N = 1, NRNU IF (I .EQ. ICKWRK(IcRNU+N-1)) IRNU = N 10 CONTINUE C c DO 100 J = 1,2 NS = 0 DO 50 N = 1, MXSP NU = ICKWRK(IcNU + (I-1)*MXSP + N - 1) K = ICKWRK(IcNK + (I-1)*MXSP + N - 1) C IF (IRNU .GT. 0) THEN RNU = RCKWRK(NcRNU + (IRNU-1)*MXSP + N - 1) IF (ABS(RNU) .GT. CKMIN) THEN IF (J .EQ. 1) THEN NU = -1 ELSE NU = 1 ENDIF ENDIF ENDIF IF (J.EQ.1.AND.NU.LT.0 .OR. J.EQ.2.AND.NU.GT.0) THEN NS = NS + 1 C IF (NS .GT. 1) THEN LT = ILASCH(ISTR) IF (LT+1 .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) = '+' ENDIF IF (IRNU .GT. 0) THEN CALL CKR2CH (ABS(RNU), IDUM, L, IERR) ELSE CALL CKI2CH (IABS(NU), IDUM, L, IERR) ENDIF IF (IERR) THEN KERR = .TRUE. WRITE (LOUT,*) ' Syntax error in CKSYMR...' ISTR = ' ' RETURN ENDIF IF (IABS(NU).NE.1 .OR. 1 (IRNU.GT.0 .AND. ABS(RNU).NE.1.0)) THEN LT = ILASCH(ISTR) IF (LT+L .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) = IDUM ENDIF LK = ILASCH(CCKWRK(IcKK+K-1)) LT = ILASCH(ISTR) IF (LT+LK .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) = CCKWRK(IcKK+K-1)(:LK) ENDIF 50 CONTINUE C DO 60 N = 1, NFAL IF (ICKWRK(IcFL+N-1) .EQ. I) THEN LT = ILASCH(ISTR) IF (ICKWRK(IcKF+N-1) .EQ. 0) THEN IF (LT+4 .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) = '(+M)' ELSE IDUM = ' ' IDUM = CCKWRK (IcKK + ICKWRK(IcKF+N-1) - 1) LK = ILASCH(IDUM) IF (LT+LK+3 .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) ='(+'//IDUM(:LK)//')' ENDIF ENDIF 60 CONTINUE C DO 70 N = 1, NTHB IF (ICKWRK(IcTB+N-1).EQ.I .AND. INDEX(ISTR,'(+M)').LE.0) 1 THEN LT = ILASCH(ISTR) IF (LT+2 .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF ISTR(LT+1:) = '+M' ENDIF 70 CONTINUE C DO 80 N = 1, NWL IF (ICKWRK(IcWL+N-1) .EQ. I) THEN W = RCKWRK(NcWL+N-1) LT = ILASCH(ISTR) IF (LT+3 .GT. ILEN) THEN KERR = .TRUE. ISTR = ' ' WRITE (LOUT, 500) RETURN ENDIF IF (J.EQ.1.AND.W.LT.0.0 .OR. J.EQ.2.AND.W.GT.0.0) 1 ISTR(LT+1:) = '+HV' ENDIF 80 CONTINUE C c IF (J.EQ.1) THEN c LT = ILASCH(ISTR) c IF (ICKWRK(IcNS+I-1) .LT. 0) THEN c IF (LT+2 .GT. ILEN) THEN c KERR = .TRUE. c ISTR = ' ' c WRITE (LOUT, 500) c RETURN c ENDIF c ISTR(LT+1:) = '=>' c ELSE c IF (LT+3 .GT. ILEN) THEN c KERR = .TRUE. c ISTR = ' ' c WRITE (LOUT, 500) c RETURN c ENDIF c ISTR(LT+1:) = '<=>' c ENDIF c ENDIF c 100 CONTINUE LT = ILASCH(ISTR) C 500 FORMAT (' Error in CKSYMR...character string length too small') RETURN END C C----------------------------------------------------------------------C C SUBROUTINE CKDUP (I, MAXSP, NS, NR, NU, NUNK, IFAL, KFAL, 1 IW, RW, ISAME) C C Checks reaction I against the (I-1) reactions for duplication C----------------------------------------------------------------------C C*****precision > double IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) C*****END precision > double C*****precision > single C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N) C*****END precision > single C DIMENSION NS(*), NR(*), NU(MAXSP,*), NUNK(MAXSP,*), IFAL(*), 1 KFAL(*), IW(*), RW(*), ithb(MAXSP), irev(NII) include 'ckstrt.h' C ISAME = 0 NRI = NR(I) NPI = ABS(NS(I)) - NR(I) C DO 500 J = 1, I-1 C NRJ = NR(J) NPJ = ABS(NS(J)) - NR(J) C IF (NRJ.EQ.NRI .AND. NPJ.EQ.NPI) THEN C NSAME = 0 DO 20 N = 1, MAXSP KI = NUNK(N,I) NI = NU(N,I) C DO 15 L = 1, MAXSP KJ = NUNK(L,J) NJ = NU(L,J) IF (NJ.NE.0 .AND. KJ.EQ.KI .AND. NJ.EQ.NI) 1 NSAME = NSAME + 1 15 CONTINUE 20 CONTINUE C IF (NSAME .EQ. ABS(NS(J))) THEN C C same products, reactants, coefficients, check fall-off C third body C c zsigy c This loop would correctly distinguish two differently parametrized c fall-off reactions. c c if (NFAL.GT.0) then c do 999 l=1, NFAL c IF (IFAL(l).EQ.I) THEN c DO 22 N = 1, NFAL c IF (J.EQ.IFAL(N) .AND. KFAL(N).EQ.KFAL(l)) THEN c ISAME = J c RETURN c ENDIF c 22 CONTINUE c ENDIF c 999 continue c endif c c The original loop from ckinterp: c zsigy IF (NFAL.GT.0 .AND. IFAL(NFAL).EQ.I) THEN DO 22 N = 1, NFAL-1 IF (J.EQ.IFAL(N) .AND. KFAL(N).EQ.KFAL(NFAL)) THEN ISAME = J RETURN ENDIF 22 CONTINUE RETURN ENDIF C c zsigy checking the third body call ckitr(iw, rw, ithb, irev) if (ithb(i).eq. -1 .and. ithb(j).eq. -1) isame=j if (ithb(i).ge. 0 .and. ithb(j).ge. 0) isame=j c zsigy c ISAME = J RETURN ENDIF ENDIF C IF (NPI.EQ.NRJ .AND. NPJ.EQ.NRI) THEN C NSAME = 0 DO 30 N = 1, MAXSP KI = NUNK(N,I) NI = NU(N,I) C DO 25 L = 1, MAXSP KJ = NUNK(L,J) NJ = NU(L,J) IF (NJ.NE.0 .AND. KJ.EQ.KI .AND. -NJ.EQ.NI) 1 NSAME = NSAME + 1 25 CONTINUE 30 CONTINUE C IF (NSAME.EQ.ABS(NS(J)) .AND. 1 (NS(J).GT.0 .OR. NS(I).GT.0)) THEN C C same products as J reactants, and vice-versa C IF (NFAL.GT.0 .AND. IFAL(NFAL).EQ.I) THEN DO 32 N = 1, NFAL-1 IF (J.EQ.IFAL(N) .AND. KFAL(N).EQ.KFAL(NFAL)) THEN ISAME = J RETURN ENDIF 32 CONTINUE RETURN ENDIF C c zsigy checking the third body call ckitr(iw, rw, ithb, irev) if (ithb(i).eq. -1 .and. ithb(j).eq. -1) isame=j if (ithb(i).ge. 0 .and. ithb(j).ge. 0) isame=j c zsigy c ISAME = J RETURN ENDIF ENDIF C 500 CONTINUE RETURN END c c--------------------------------------------------------------------- c CHARACTER*(*) FUNCTION UPCAS(ISTR, ILEN) CHARACTER ISTR*(*), LCASE(26)*1, UCASE(26)*1 DATA LCASE /'a','b','c','d','e','f','g','h','i','j','k','l','m', 1 'n','o','p','q','r','s','t','u','v','w','x','y','z'/, 2 UCASE /'A','B','C','D','E','F','G','H','I','J','K','L','M', 3 'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/ C UPCAS = ' ' UPCAS = ISTR(:ILEN) JJ = MIN (LEN(UPCAS), LEN(ISTR), ILEN) DO 10 J = 1, JJ DO 10 N = 1,26 IF (ISTR(J:J) .EQ. LCASE(N)) UPCAS(J:J) = UCASE(N) 10 CONTINUE RETURN END c c--------------------------------------------------------------------- c subroutine control(lctr,lout,lnul,kk,kname,L16,iformat,ltherm, 1 lrev,metrev,nkill,mxkill,ikill,iunitA,iunitE, 2 ithinfo, mxmod, nHmod, iHmod, rHmod, nSmod, iSmod, rSmod, 3 Tlow,Tmid,Thigh) c c Keyword found: output c IRREV {n} lrev= .false. (default: .true.) c metrev= n (default: n=1 ) c c THERMO ltherm= .true. (default: .false.) c c KILL nkill= number of species to be eliminated c ikill(.)= list of species to be eliminated c c FORMAT1 iform=1 (default) c FORMAT2 iform=2 c c MOLES iunitA=1 (default) c MOLECULES iunitA=2 c c KELVINS iunitE=1 (default) c CAL/MOLE iunitE=2 c KCAL/MOLE iunitE=3 c JOULES/MOLE iunitE=4 c KJOULES/MOLE iunitE=5 c c ithinfo=0 (default) c THINFO_KJ ithinfo=1 c THINFO_KCAL ithinfo=2 c c c Thermodynamic data modifying commands: c c nHmod= number of species with modified dHf c NEWH_KJ | iHmod(.)= species_number c NEWH_KCAL | dHnew(.)= new_data c c nSmod= number of species with modified S c NEWS_J/K | iSmod(.)= species_number c NEWS_CAL/K | dSnew(.)= new_data c implicit double precision (a-h, o-z), integer (i-n) parameter (nc=19) dimension nray(nc), ikill(mxkill), iHmod(mxmod), iSmod(mxmod) dimension rHmod(mxmod), rSmod(mxmod) character*(*) kname(kk) character*80 line, upcas character*16 kray(nc) logical lrev, ltherm, kerr, L16 c c------------------------------------------- c c Commands c kray( 1)= 'IRREV' kray( 2)= 'FORMAT1' kray( 3)= 'FORMAT2' kray( 4)= 'MOLES' kray( 5)= 'MOLECULES' kray( 6)= 'KELVINS' kray( 7)= 'CAL/MOLE' kray( 8)= 'KCAL/MOLE' kray( 9)= 'JOULES/MOLE' kray(10)= 'KJOULES/MOLE' kray(11)= 'KILL' kray(12)= 'THERMO' kray(13)= 'THINFO_KJ' kray(14)= 'THINFO_KCAL' kray(15)= 'NEWH_KJ' kray(16)= 'NEWH_KCAL' kray(17)= 'NEWS_J/K' kray(18)= 'NEWS_CAL/K' kray(19)= 'END' c c defaults c lrev= .true. ltherm= .false. metrev= 1 nkill= 0 nkill1= 0 iformat= 1 iunitA= 1 iunitE= 1 ithinfo= 0 nHmod= 0 nSmod= 0 c write(lout,200) 200 format(//5x,'The control file: '/) c c===================================================== c 100 continue line = ' ' read (lctr,'(a)',end=99) line write(lout,'(a)') line ilen = ipplen(line) c c blank lines and comment lines c if (ilen .eq. 0) go to 100 line= upcas(line,ilen) c c commands c call ckcray(line(:ilen),nc,kray,lnul,nc,nray,nf,kerr) if (nf.eq.0) then write(lout,201) 201 format(/' Sorry, I could not understand this line.'/) goto 100 endif if (nf.gt.1) then write(lout,202) 202 format(1x,'*** Please use only one command on each line.', 1 1x,'This line is ignored. ***') goto 100 endif icid= nray(1) goto (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,99) icid goto 100 c c IRREV c 1 lrev= .false. n= 1 call cksnum(line(:ilen),n,lnul,kray,nc,knum,nval,rval,kerr) if (nval.eq.0) goto 100 if (nval.gt.1) then write(lout,224) 224 format(' The numbers are ignored.'/) goto 100 endif mr= idint(rval) if (mr.eq.1 .or. mr.eq.2 .or. mr.eq.3) then metrev= mr else write(lout,203) mr 203 format (' Sorry, there is no reversal method No',i5/) endif goto 100 c c FORMAT1 c 2 iformat= 1 goto 100 c c FORMAT2 c 3 iformat= 2 goto 100 c c MOLES c 4 iunitA= 1 goto 100 c c MOLECULES c 5 iunitA= 2 goto 100 c c KELVINS c 6 iunitE= 1 goto 100 c c CAL/MOLE c 7 iunitE= 2 goto 100 c c KCAL/MOLE c 8 iunitE= 3 goto 100 c c JOULES/MOLE c 9 iunitE= 4 goto 100 c c KJOULES/MOLE c 10 iunitE= 5 goto 100 c c KILL c 11 continue c call ckcray(line,kk,kname,lnul,kk,ikill(nkill+1),nkill1,kerr) if (nkill1.eq.0) then write(lout,204) 204 format(1x,'*** I could find a single species'' name. ', 1 1x,'This line is ignored ***') endif nkill= nkill+nkill1 goto 100 c c THERMO c 12 continue ltherm= .true. goto 100 c c THINFO_KJ c 13 continue ithinfo= 1 ltherm= .true. goto 100 c c THINFO_KCAL c 14 continue ithinfo= 2 ltherm= .true. goto 100 c c NEWH_KJ c 15 continue nHmod= nHmod+1 line= line(ifirch(line)+8:) call cksnum(line, 1, lnul, kname, kk, 1 iHmod(nHmod), nval, rHmod(nHmod), kerr) ltherm= .true. goto 100 c c NEWH_KCAL c 16 continue nHmod= nHmod+1 line= line(ifirch(line)+10:) call cksnum(line, 1, lnul, kname, kk, 1 iHmod(nHmod), nval, rHmod(nHmod), kerr) ltherm= .true. c c conversion from kcal->kJ c rHmod(nHmod)= rHmod(nHmod) * 4.184 goto 100 c c NEWS_J/K c 17 continue nSmod= nSmod+1 line= line(ifirch(line)+9:) call cksnum(line, 1, lnul, kname, kk, 1 iSmod(nSmod), nval, rSmod(nSmod), kerr) ltherm= .true. goto 100 c c NEWS_CAL/K c 18 continue nSmod= nSmod+1 line= line(ifirch(line)+11:) call cksnum(line, 1, lnul, kname, kk, 1 iSmod(nSmod), nval, rSmod(nSmod), kerr) ltherm= .true. c c conversion from cal/K->J/K c rSmod(nSmod)= rSmod(nSmod) * 4.184 goto 100 c c---------------------------------------------------------------- c c the accepted commands c 99 continue write(lout,205) 205 format(//' Features of the transformed mechanism: '//) c c REVERSIBILITY c if (lrev) then write(lout,206) 206 format(' Reversible reactions will be kept reversible'/) else write(lout,207) 207 format(' All reversible reactions will be transformed to pairs of' 1 ,' irreversible reactions'/) c goto (25,26,27) metrev goto 28 25 write(lout,208) Tlow, Tmid, Thigh 208 format(' Reverse Arrhenius parameters A, n, and E will be fitted'/ 1 ' at the following 3 temperatures:'/ 2 f7.1,' K, ',f6.1,' K, and ',f6.1,' K'/) goto 28 26 write(lout,209) Tlow, Thigh 209 format( 1 ' Forward Arrhenius parameter n will be preserved as reverse n'/ 2 ' Reverse Arrhenius parameters A and E will be fitted'/ 3 ' at the following 2 temperatures:'/ 4 f7.1,' K and ',f6.1,' K'/) goto 28 27 write(lout,210) Tlow, Thigh 210 format( 1 ' Forward Arrhenius parameter A will be preserved as reverse A'/ 2 ' Reverse Arrhenius parameters n and E will be fitted'/ 3 ' at the following 2 temperatures:'/ 4 f7.1,' K and ',f6.1,' K'/) 28 continue endif c c THINFO c goto (50,51,52) (ithinfo+1) goto 53 50 write(lout,250) 250 format(' Thermodynamic summary info will not be printed'/) goto 53 51 write(lout,251) 251 format(' Thermodynamic summary info will be printed', 1 ' in Joule units'/) goto 53 52 write(lout,252) 252 format(' Thermodynamic summary info will be printed', 1 ' in calorie units'/) 53 continue c c NEW H c if (nHmod.eq.0) then write(lout,260) 260 format(' No enthalpy value is modified') goto 69 endif c write(lout,261) 261 format(//' Modified enthalpies; new values:'/ 1 ' species dH (kJ) dH (kcal)') do 60 i=1,nHmod write(lout,262) kname(iHmod(i)), rHmod(i), rHmod(i)/4.184 262 format(1x,A16,f8.3,3x,f8.3) 60 continue 69 continue c c NEW S c if (nSmod.eq.0) then write(lout,270) 270 format(/' No entropy value is modified') goto 79 endif c write(lout,271) 271 format(//' Modified entropies'/ 1 ' species S (J/K) S (cal/K)') do 70 i=1,nSmod write(lout,262) kname(iSmod(i)), rSmod(i), rSmod(i)/4.184 70 continue 79 continue c c THERMO c if (ltherm) then write(lout,225) else write(lout,226) endif 225 format 1 (/' Thermodynamic data will be printed in the mechanism file'/) 226 format 1 (/' Thermodynamic data will not appear in the mechanism file'/) c c FORMAT c goto (21,22) iformat goto 23 21 write(lout,211) 211 format 1 (' Reactions in the mechanism file will look like this:'/ 2 ' CH4 + OH => CH3 + H2O 1.570E07 1.8300 11.64'/) c goto 23 22 write(lout,212) 212 format 1 (' Reactions in the mechanism file will look like this:'/ 2 ' CH4 + OH => CH3 + H2O 1.570E07 1.830 11.64'/) c 23 continue c c UNITS OF A c goto (31,32) iunitA goto 33 31 write(lout,213) 213 format(' Units for A will be moles-cm-sec-K'/) goto 33 32 write(lout,214) 214 format(' Units for A will be molecules-cm-sec-K'/) 33 continue c c UNITS FOR E c goto (41,42,43,44,45) iunitE goto 46 41 write(lout,215) 215 format(' Units for E will be Kelvin'/) goto 46 42 write(lout,216) 216 format(' Units for E will be cal/mole'/) goto 46 43 write(lout,217) 217 format(' Units for E will be kcal/mole'/) goto 46 44 write(lout,218) 218 format(' Units for E will be Joules/mole'/) goto 46 45 write(lout,219) 219 format(' Units for E will be kJoules/mole'/) 46 continue c c REACTION ELIMINATION c if(nkill.eq.0) then write(lout,220) 220 format(' No reaction will be eliminated'//) return endif write(lout,221) 221 format(' All consuming and producing reactions of the '/ 1 ' following species will be eliminated: ') if (L16) then write(lout,222) (kname(ikill(k)),k=1,nkill) 222 format(5A16) else write(lout,223) (kname(ikill(k)),k=1,nkill) 223 format(10a8) endif c return end c c-------------------------------------------------------------------------- c subroutine rwrite 1 (lout,lmech,iformat,lrev,lirrev,istrl,istrr,A,Tn,E,iferr,NoREA) implicit double precision (a-h, o-z), integer (i-n) logical lrev,lirrev character*(*) istrl,istrr character*80 line iferr= 0 ill= ilasch(ISTRL) ilr= ilasch(ISTRR) if (iformat.eq.1) then c c format 1 output c if (lrev .and. .not.lirrev ) then line= ISTRL(:ill) // ' = ' // ISTRR else line= ISTRL(:ill) // ' => ' // ISTRR endif c write(lmech,306) line(:49),A,Tn,E 306 format(a49,1x,1pe10.3,1x,0pf8.4,1x,f10.2) else c c format 2 output c if(ill.gt.23 .or. ilr.gt.23) then iformat=1 write(lout,210) NoREA 210 format(5x,' The species names in reaction', i6,' are too long'/ 1 5x,' for format 2. Switch to format 1'//) c iferr= 1 return else if (lrev .and. .not.lirrev ) then write(lmech,307) ISTRL,ISTRR,A,Tn,E 307 format(a23,' = ',a23,1x,1pe10.3,1x,0pf7.3,1x,f10.2) else write(lmech,308) ISTRL,ISTRR,A,Tn,E 308 format (a23,' => ',a23,1x,1pe10.3,1x,0pf7.3,1x,f10.2) endif endif endif return end c c-------------------------------------------------------------------- c subroutine backA(metrev,i,iw,rw,SMH,T1,T2,T3, 1 A,Tn,E,Ar,Tnr,Er,rvers,kk) implicit double precision (a-h, o-z), integer (i-n) COMMON /MACH/ SMALL,BIG,EXPARG dimension iw(1),rw(1),SMH(KK) c c include 'ckstrt.h' NPAR4=NPAR+1 c Ru= rw(NcRU) Patm= rw(NcPA) c c the forward k c fk1 = A * dexp(Tn * dlog(T1) - E/T1) fk2 = A * dexp(Tn * dlog(T2) - E/T2) fk3 = A * dexp(Tn * dlog(T3) - E/T3) c c the equilibrium constant at T1 c CALL CKSMH (T1, iw, rw, SMH) SUMSMH = 0.0 nusumk = 0 DO 1 N = 1, MXSP NUNKNI= iw(IcNK+(i-1)*MXSP+n-1) nuni= iw(IcNU+(i-1)*MXSP+n-1) IF (NUNKNI.NE.0) SUMSMH=SUMSMH+nuni*SMH(NUNKNI) nusumk= nusumk+nuni 1 CONTINUE eqK = dexp(MIN(SUMSMH,EXPARG)) c c the reverse k at T1 c PFAC = Patm / (Ru*T1) eqK = eqK * PFAC**nusumk rk1 = fk1 / MAX(eqK,SMALL) c c the equilibrium constant at T2 c CALL CKSMH (T2, iw, rw, SMH) SUMSMH = 0.0 nusumk = 0 DO 2 N = 1, MXSP NUNKNI= iw(IcNK+(i-1)*MXSP+n-1) nuni= iw(IcNU+(i-1)*MXSP+n-1) IF (NUNKNI.NE.0) SUMSMH=SUMSMH+nuni*SMH(NUNKNI) nusumk= nusumk+nuni 2 CONTINUE eqK = dexp(MIN(SUMSMH,EXPARG)) c c the reverse k at T2 c PFAC = Patm / (Ru*T2) eqK = eqK * PFAC**nusumk rk2 = fk2 / MAX(eqK,SMALL) c c the equilibrium constant at T3 c CALL CKSMH (T3, iw, rw, SMH) SUMSMH = 0.0 nusumk = 0 DO 3 N = 1, MXSP NUNKNI= iw(IcNK+(i-1)*MXSP+n-1) nuni= iw(IcNU+(i-1)*MXSP+n-1) IF (NUNKNI.NE.0) SUMSMH=SUMSMH+nuni*SMH(NUNKNI) nusumk= nusumk+nuni 3 CONTINUE eqK = dexp(MIN(SUMSMH,EXPARG)) c c the reverse k at T3 c PFAC = Patm / (Ru*T3) eqK = eqK * PFAC**nusumk rk3 = fk3 / MAX(eqK,SMALL) c c Method 1: Fitting A,n, and E to reverse k values c at three temperatures T1, T2, and T3 c c Method 2: Retain forward Tn and c calculate A and E knowing reverse k at T1 and T3 c c Method 3: Retain forward A and c calculate Tn and E knowing reverse k at T1 and T3 c goto (10,20,30) metrev 10 continue f1= dlog(rk1) f2= dlog(rk2) f3= dlog(rk3) a1= dlog(T1) a2= dlog(T2) a3= dlog(T3) c den= -T2*a1*T1+T2*a3*T3+T3*a1*T1+T1*a2*T2-T1*a3*T3-T3*a2*T2 Tnr= (-T2*f1*T1-T1*f3*T3+T2*f3*T3-T3*f2*T2+T3*f1*T1+T1*f2*T2)/den Aln= -(-T2*a2*f1*T1+T2*a2*f3*T3+T2*f2*a1*T1-T2*f2*a3*T3 1 +a3*T3*f1*T1-f3*T3*a1*T1)/den Ar= dexp(Aln) Er= -T3*T2*T1*(-f3*a1+a2*f3-a2*f1+f2*a1+a3*f1-a3*f2)/den return 20 continue Tnr= Tn f1= dlog(rk1) f3= dlog(rk3) a1= dlog(T1) a3= dlog(T3) Er= (f1-f3-Tn*a1+Tn*a3)/( 1.D0/T3-1.D0/T1) Ar= dexp( f3-Tn*a3+Er/T3 ) return 30 continue Ar= A Aln= dlog(A) f1= dlog(rk1) f3= dlog(rk3) a1= dlog(T1) a3= dlog(T3) Tnr= (T1*f1-T2*f3-T1*Aln+T2*Aln)/(T1*a1-T2*a3) Er= T3*(Aln+Tnr*a3-f3) return end c c-------------------------------------------------------------------- c subroutine backLT (i,iw,rw,SMH,T1,T3, 1 A,Tn,E,TL1,TL2,Ar,Tnr,Er,TL1r,TL2r,rvers,kk) implicit double precision (a-h, o-z), integer (i-n) COMMON /MACH/ SMALL,BIG,EXPARG dimension iw(1),rw(1),SMH(KK) c c include 'ckstrt.h' NPAR4=NPAR+1 c Ru= rw(NcRU) Patm= rw(NcPA) c c the forward k c fk1 = A * dexp(Tn * dlog(T1) - E/T1) 1 * dexp( TL1*T1**(-1./3.) + TL2*T1**(-2./3.) ) fk3 = A * dexp(Tn * dlog(T3) - E/T3) 1 * dexp( TL1*T3**(-1./3.) + TL2*T3**(-2./3.) ) c c The equilibrium constant at T1 c CALL CKSMH (T1, iw, rw, SMH) SUMSMH = 0.0 nusumk = 0 DO 1 N = 1, MXSP NUNKNI= iw(IcNK+(i-1)*MXSP+n-1) nuni= iw(IcNU+(i-1)*MXSP+n-1) IF (NUNKNI.NE.0) SUMSMH=SUMSMH+nuni*SMH(NUNKNI) nusumk= nusumk+nuni 1 CONTINUE eqK = dexp(MIN(SUMSMH,EXPARG)) c c the reverse k at T1 c PFAC = Patm / (Ru*T1) eqK = eqK * PFAC**nusumk rk1 = fk1 / MAX(eqK,SMALL) c c Removal of Arrhenius' contribution c rk1= rk1/( Ar * dexp(Tnr * dlog(T1) - Er/T1) ) c c The equilibrium constant at T3 c CALL CKSMH (T3, iw, rw, SMH) SUMSMH = 0.0 nusumk = 0 DO 3 N = 1, MXSP NUNKNI= iw(IcNK+(i-1)*MXSP+n-1) nuni= iw(IcNU+(i-1)*MXSP+n-1) IF (NUNKNI.NE.0) SUMSMH=SUMSMH+nuni*SMH(NUNKNI) nusumk= nusumk+nuni 3 CONTINUE eqK = dexp(MIN(SUMSMH,EXPARG)) c c The reverse k at T3 c PFAC = Patm / (Ru*T3) eqK = eqK * PFAC**nusumk rk3 = fk3 / MAX(eqK,SMALL) c c Removal of Arrhenius' contribution c rk3= rk3/( Ar * dexp(Tnr * dlog(T3) - Er/T3) ) c f1= dlog(rk1) f3= dlog(rk3) a1= 1./T1**(1./3.) a3= 1./T3**(1./3.) b1= 1./T1**(2./3.) b3= 1./T3**(2./3.) TL2r= (f1/a1-f3/a3)/(b1/a1-b3/a3) TL1r= (f3-TL2r*b3)/a3 return end c c-------------------------------------------------------------------- c subroutine convert(i,MXSP,NU,ico,ifltyp,l3rd, 1 iunitA,iunitE,A,E,Ac,Ec,rvers) implicit double precision (a-h, o-z), integer (i-n) dimension NU(MXSP,1) logical l3rd c c Conversion of A units depends on the following status flags: c c ico 1 Forward Arrhenius data c 2 Backward Arrhenius data c 3 LOW pressure limit - forward c 4 LOW pressure limit - backward c c ifltyp =0 not fall-off reaction c >0 fall-off reaction c c l3rd .true. Includes 3rd body c c if (rvers.lt.3.7999) then if (iunitE.eq.1) EFAC= 1.000 if (iunitE.eq.2) EFAC= 1.987 if (iunitE.eq.3) EFAC= 1.987 / 1000. if (iunitE.eq.4) EFAC= 8.314 if (iunitE.eq.5) EFAC= 8.314 / 1000. else if (iunitE.eq.1) EFAC= 1.000 if (iunitE.eq.2) EFAC= 8.314510D0 / 4.184d0 if (iunitE.eq.3) EFAC= 8.314510D-3 / 4.184d0 if (iunitE.eq.4) EFAC= 8.314510D0 if (iunitE.eq.5) EFAC= 8.314510D-3 endif Ec= E * EFAC c Ac= A if (iunitA.eq.1) return c NSTOR = 0 NSTOP = 0 do 5 n = 1, MXSP if (NU(n,i) .lt. 0) then c sum of stoichiometric coefficients of reactants NSTOR = NSTOR + abs(NU(n,i)) elseif (NU(n,i) .gt. 0) then c sum of stoichiometric coefficients of products NSTOP = NSTOP + NU(n,i) endif 5 continue c AVAG = 6.023E23 c goto (1,2,3,4) ico c c Forward Arrhenius data c 1 if (NSTOR .le. 0) goto 9 if (l3rd .and. ifltyp.eq.0) then c 3rd body and not fall-off Ac = A / AVAG**NSTOR else c no 3rd body or fall-off Ac = A / AVAG**(NSTOR-1) endif goto 9 c c Backward Arrhenius data c 2 if (NSTOP .le. 0) goto 9 if (l3rd .and. ifltyp.eq.0) then c 3rd body and not fall-off Ac = A / AVAG**NSTOP else c no 3rd body or fall-off Ac = A / AVAG**(NSTOP-1) endif goto 9 c c LOW Arrhenius data - forward c 3 if (NSTOR .le. 0) goto 9 Ac = A / AVAG**NSTOR goto 9 c c LOW Arrhenius data - backward c 4 if (NSTOP .le. 0) goto 9 Ac = A / AVAG**NSTOP c 9 continue return end c c-------------------------------------------------------------------- c subroutine r2ch(r,ch,lch,kerr) implicit double precision (a-h, o-z), integer (i-n) character*5 ch, chi logical kerr c c conversion of real numbers (0<=r<99.995) to a c max. 5-character string of 'nn.nn' format c if (r.gt.99.995) then kerr=.true. return endif ir= r*100. call cki2ch (ir, chi, l, kerr) if (kerr) return if (l.eq.1) ch='0.0'//chi(:l) if (l.eq.2) ch='0.' //chi(:l) if (l.gt.2) ch=chi(:(l-2))//'.'//chi( (l-1) : l ) lch= ilasch(ch) if (ch((lch-1):lch).eq.'00') ch=ch(1:(lch-3)) if (ch(lch:lch) .eq.'0' ) ch=ch(1:(lch-1)) lch= ilasch(ch) return end c c 1 2 3 4 5 6 7 c23456789012345678901234567890123456789012345678901234567890123456789012