subroutine chem_init c----------------------------------------------------------------------------- IMPLICIT NONE include 'parameter.par' include 'chem.com' integer mdim,kdim,mkdim,lsym,npidim,npc,npcp2,maxtp,ntr, 1 nktdim,nlar,nsidim,ntidim,nlidim,nfidim,ntdim, 2 nidim,lin,lout,lthrm,linc,ckmin,maxord,noidim parameter (mdim=10, kdim=nsmx, mkdim=mdim*kdim, lsym=16, 1 npidim=idim*npar, npc=5, npcp2=npc+2, maxtp=3, 2 ntr=maxtp-1, nktdim=ntr*nPcp2*kdim, 3 nlar=2, nsidim=maxsp*idim, ntidim=maxtb*idim, 4 nlidim=nlar*idim, nfidim=nfar*idim, 5 ntdim=kdim*maxtp, nidim=11*idim, lin=15, lout=16, 6 lthrm=17, linc=25, ckmin=1.0d-3, maxord=kdim, 7 noidim=maxord*idim) integer njar,nf1r,njidim,nf1idim parameter (njar=9, nf1r=4, njidim=njar*idim, nf1idim=nf1r*idim) character kname(kdim)*(lsym) character ename(mdim)*(lsym), 1 key(5)*4, vers*(lsym), prec*(lsym) real*8 awt(mdim), wtm(kdim), 1 a(npcp2,ntr,kdim), t(maxtp,kdim),plan(nlar,idim), 2 rlan(nlar,idim),wl(idim),rnu(maxsp,idim), 3 rord(maxord,idim),pjan(njar,idim),pft1(nf1r,idim), 4 pexc(idim) integer kncf(mdim,kdim),kphse(kdim),kchrg(kdim),nt(kdim), 1 idup(idim),ilan(idim),irlt(idim),iwl(idim), 2 irnu(idim),iord(idim),kord(maxord,idim) logical kerr, thermo, ithrm(kdim) integer ieim(idim), itdep(idim), ijan(idim),ift1(idim), 1 iexc(idim) character*80 fmt integer neim,njan,nft1,nexc data neim,njan,nft1,nexc/4*0/, ieim/idim*0/, itdep/idim*0/, 1 ijan/idim*0/, ift1/idim*0/, pjan/njidim*0.0/, 2 pft1/nf1idim*0.0/, pexc/idim*0.0/ integer itask,nchrg,nlan,nrlt,nwl, 1 nrnu,nord c data key/'elem','spec','ther','reac','end'/, kerr/.false./, c 1 itask,nchrg,mm,kk,ii,nlan,nfal,nthb,nrev,nrlt,nwl, c * nrnu,nord/13*0/, c 2 ename,awt/mdim*' ',mdim*0.0/, thermo/.true./, c 3 t/ntdim*-1.0/, kname,wtm,nt,kphse,kchrg,ithrm c 4 /kdim*' ', kdim*0.0, kdim*3, kdim*0, kdim*0, kdim*.false./, c 5 wl,ifop,ntbs,idup /idim*0.0, idim*-1, idim*0, idim*0/, c 6 nspec,nreac,irev,ilan,irlt,iwl,ifal,kfal,ithb,irnu,iord c 7 /nidim*0/ data nunk,nu/nsidim*0, nsidim*0/, nktb,aik/ntidim*0,ntidim*-1.0/ data rnu/nsidim*0.0/, kord/noidim*0/, rord/noidim*0.0/ data par,rpar/npidim*0.0, npidim*0.0/ data plan,rlan/nlidim*0.0, nlidim*0.0/ data fpar/nfidim*0.0/, kncf/mkdim*0.0/, a/nktdim*0.0/ integer lenick, lenrck, lencck, mm, kk, ii, maxsp2,maxtb2, 1 maxtp2, npc2, npar2, nlar2, nfar2, njar2, 2 nf1r2, maxord2, ckmin2 integer m,k,l,i,n,nr,ip3,kthb,kfl,nspecmx,nk,ir double precision dflt c----------------------------------------------------------------------------- open(unit=10,file='chem.bin',status='old',form='unformatted') c----------------------------------------------------------------------------- read (10) vers, prec, kerr read (10) lenick, lenrck, lencck, mm, kk, ii, maxsp2, 1 maxtb2, maxtp2, npc2, npar2, nlar2, nfar2, nrev, nfal, 2 nthb, nlan, nrlt, nwl, nchrg, neim, njar2, njan, 3 nf1r2, nft1, nexc, nrnu, nord, maxord2, ckmin2 read (10) (ename(m), awt(m), m = 1, mm) read (10) (kname(k), (kncf(m,k),m=1,mm), kphse(k), 1 kchrg(k), wtm(k), nt(k), (t(l,k),l=1,maxtp), 2 ((a(m,l,k), m=1,npcp2), l=1,ntr), k = 1, kk) c----------------------------------------------------------------------------- if (ii .gt. 0) then read (10) (nspec(i), nreac(i), (par(n,i), n = 1, npar), 1 (nu(m,i), nunk(m,i), m = 1, maxsp), i = 1, ii) if (nrev .gt. 0) read (10) 1 (irev(n),(rpar(l,n),l=1,npar),n=1,nrev) if (nfal .gt. 0) read (10) 1 (ifal(n),ifop(n),kfal(n),(fpar(l,n),l=1,nfar), n = 1, nfal) if (nthb .gt. 0) read (10) 1 (ithb(n),ntbs(n),(nktb(m,n),aik(m,n),m=1,maxtb),n=1,nthb) if (nlan .gt. 0) read (10) 1 (ilan(n), (plan(l,n), l = 1, nlar), n = 1, nlan) if (nrlt .gt. 0) read (10) 1 (irlt(n), (rlan(l,n), l = 1, nlar), n=1,nrlt) if (nwl .gt. 0) read (10) (iwl(n), wl(n), n = 1, nwl) if (neim .gt. 0) read (10) (ieim(n),itdep(n),n=1,neim) if (njan .gt. 0) read (10) 1 (ijan(n), (pjan(l,n), l = 1, njar), n = 1, njan) if (nft1 .gt. 0) read (10) 1 (ift1(n), (pft1(l,n), l = 1, nf1r), n = 1, nft1) if (nexc .gt. 0) read (10)(iexc(n), pexc(n), n=1, nexc) if (nrnu .gt. 0) read (10) 1 (irnu(n), (rnu(m,n), m = 1, maxsp), n = 1, nrnu) if (nord .gt. 0) read (10) 1 (iord(n), (kord(l,n), rord(l,n), l=1, maxord), n=1,nord) else write(6,*)' gckwk: input file has no reaction information ' endif c----------------------------------------------------------------------------- close(10) c----------------------------------------------------------------------------- #ifdef VERBOSE c----------------------------------------------------------------------------- c ... write out Arrhenius parameters open(unit=20,file='wdat.out',status='unknown') write(20,'(''# i, B, exponent, E/R0 (K) '')') do i=1,ii write(20,'(i5,3e22.14)')i,(par(k,i),k=1,3) end do close(20) c----------------------------------------------------------------------------- c ... write out reverse Arrhenius parameters for reactions with c ... specified reverse Arrhenius parameters if(nrev.ne.0)then open(unit=20,file='wrevdat.out',status='unknown') write(20,'(''# i, irev, B, exponent, E/R0 (K) '')') do i=1,nrev write(20,'(2i5,3e22.14)')i,irev(i),(rpar(k,i),k=1,3) end do close(20) endif c----------------------------------------------------------------------------- c ... write out species_in_reaction data open(unit=20,file='spinr.out',status='unknown') write(20,'(''# number of species: kk = '',i5)')kk do k=1,kk write(20,'(''# '',i5,5x,a16)')k,kname(k) end do write(20,'(''# max allowed spec/reacton: maxsp = '',i5)')maxsp write(20,'(''# nspec(i) = number of species in reaction i'')') write(20,'(''# nreac(i) = number of reactants in reacn i'')') write(20,'(''# nu(k,i) = stoich coef species k in reacn i'')') write(20,'(''# nunk(k,i)= index (1...kk) specie k reacn i'')') write(20,'(''# i,nspec,nreac,(nu(k,i),nunk(k,i),k=1,maxsp)'')') do i=1,ii write(20,*)i,nspec(i),nreac(i), + (nu(k,i),nunk(k,i),k=1,maxsp) end do close(20) c----------------------------------------------------------------------------- c ... write out third body data open(unit=20,file='thb.out',status='unknown') write(20,'(''# total number of reactions: ii='',i5)')ii write(20,'(''# num reac w/ 3rd bodies: nthb='',i5)')nthb write(20,'(''# ithb(j) (j=1..nthb) = '', + '' reaction index of j (1...ii)'')') write(20,'(''# ntbs(j) = no. enhanced species in reacn j'')') write(20,'(''# max allowed num enhanced species: maxtb='',i5)') + maxtb write(20,'(''# nktb(k,j) = species indices of enhanced species'' + ,'' k (k=1...maxtb)'')') write(20,'(''# in enhanced reaction j (1...nthb)'')') write(20,'(''# j, ithb(j), ntbs(j), (nktb(k,j),k=1,maxtb)'')') do j=1,nthb write(20,*)j,ithb(j),ntbs(j),(nktb(k,j),k=1,maxtb) end do close(20) open(unit=20,file='facthb.out',status='unknown') write(20,'(''# total number of reactions: ii='',i5)')ii write(20,'(''# num reac w/ 3rd bodies: nthb='',i5)')nthb write(20,'(''# max allowed num enhanced species: maxtb='',i5)') + maxtb write(20,'(''# aik(k,j) = Enhancement factor of enhanced '' + ,''species k (k=1...maxtb)'')') write(20,'(''# in enhanced reaction j (1...nthb)'')') write(20,'(''# j, (aik(k,j),k=1,maxtb)'')') write(fmt,'(a,i3,a)') '(i3,2x,', maxtb, 'f7.3)' do j=1,nthb write(20,fmt)j,(aik(k,j),k=1,maxtb) end do close(20) c----------------------------------------------------------------------------- #endif c----------------------------------------------------------------------------- c----------------------------------------------------------------------------- c ... fill up my arrays c------------------------------------------------------------------------------- c ... for reaction k, there are nreac(k) reactants (max 3) c ... and therefore nreac(k) products (max 3) c ... nu(1,k) ... nu(nreac(k) , k) are the negatives of the c ... stoichiometric coefficients of the reactants c ... nu(3+1,k) ... nu(3+nreac(k) , k) are the stoichiometric c ... coefficients of the products c ... nunk(i,k) is the index in ksp=1,2,...,ns vector of species c ... for i-th spcies in reaction k c ... (i=1,...,nreac(k),3+1,...,3+nreac(k)) c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... spname(k) = name of species k c ... wt(k) = dimensional molar weight of species k (kg/kg-mol) ns = kk ! number of species numreac = ii ! number of reactions nr = numreac do k=1,ns spname(k) = kname(k) wt(k) = wtm(k) end do c------------------------------------------------------------------------------- c ... nprod(k) = number of products in reaction k do k=1,numreac nprod(k) = abs(nspec(k))-nreac(k) end do c------------------------------------------------------------------------------- c ... nursum(k)=Sum_i=1,N nu_ik^{prime}, where N = total number of species c ... (if a species is not a reactant in present reaction, it c ... has 0-contribution to the sum) do k=1,numreac nursum(k) = 0 do i=1,nreac(k) nursum(k) = nursum(k) - nu(i,k) end do end do c------------------------------------------------------------------------------- c ... nupsum(k)=Sum_i=1,N nu_ik^{double-prime} do k=1,numreac nupsum(k) = 0 do i=1,nprod(k) ip3=i+3 nupsum(k) = nupsum(k) + nu(ip3,k) end do end do c------------------------------------------------------------------------------- c ... nb. a reaction such as : H + H + M = H2 + M c ... has nspec=2, nreac=1,nprod=1 c ... hence: 1-the H and H in the reactants are simply : 2H c ... 2-the third body does not appear as a 'species' c ... Also, for every species involved in the reaction k c ... even when repeated once as reactant and once as product c ... nu_ik = nu_ik^{double_prime} - nu_ik^{prime} = nu(i,k) c ... where nu(i,k) is < 0 for a reactant, > 0 for a product c------------------------------------------------------------------------------- c ... nusum(k) = Sum_i=1,N nu_ik do k=1,numreac nusum(k) = 0 do i=1,nreac(k) nusum(k) = nusum(k) + nu(i,k) end do do i=1,nprod(k) ip3=i+3 nusum(k) = nusum(k) + nu(ip3,k) end do end do c------------------------------------------------------------------------------- c ... work-saving arrays do k=1,numreac anusum(k)=dflt(nusum(k)) do i=1,maxsp anu(i,k) = dflt(nu(i,k)) end do end do c------------------------------------------------------------------------------- c ... define flag for each reaction k=1,2,...,numreac c ... jthb(k)=0 => reaction k is not a third body reaction c ... jthb(k)!=0 => reaction k is a third body reaction c ... and jthb(k) is kthb, the index of this reaction in c ... the list of 3rd body reactions do k=1,numreac jthb(k)=0 end do if(nthb.gt.0)then do kthb=1,nthb k=ithb(kthb) jthb(k)=kthb end do endif c------------------------------------------------------------------------------- c ... define flag for each reaction k=1,2,...,numreac c ... itfdep(k)=0 => forward reaction Arrhenius expression not f(T) c ... i.e. both alpha_k and E_k are 0.0 c ... itfdep(k)!=0 => forward reaction Arrhen exp is f(T) c ... i.e. either alpha_k or E_k is non-zero c ... for all reactions with specified reverse Arrhenius parameters: c ... itrdep(k)=0 => rev reaction Arrh exp not f(T) c ... itrdep(k)!=0 => rev reaction Arrh exp is f(T) c ... and, for all reactions with no-specified reverse Arrhenius c ... parameters, itrdep(k)=0, by default (not used) do k=1,numreac itfdep(k)=0 itrdep(k)=0 end do do k=1,numreac if(par(2,k).ne.0.0.or.par(3,k).ne.0.0)then itfdep(k)=1 endif end do if(nrev.ne.0)then do k=1,nrev if(rpar(2,k).ne.0.0.or.rpar(3,k).ne.0.0)then itrdep(irev(k))=1 endif end do endif c------------------------------------------------------------------------------- c ... define flag for each reaction k=1,2,...,numreac c ... jfal(k)=0 => reaction k is not a pressure fall-off reaction c ... jfal(k)!=0 => reaction k is a pressure fall-off reaction c ... and jfal(k) is kfl, the index of this reaction in c ... the list of 3rd body reactions do k=1,numreac jfal(k)=0 end do if(nfal.gt.0)then do kfl=1,nfal k=ifal(kfl) jfal(k)=kfl end do endif c------------------------------------------------------------------------------- if(numreac.gt.idim)then stop 'gckwk: numreac > idim' endif if(nrev.gt.idim)then stop 'gckwk: nrev > idim' endif nspecmx = 0 do k=1,numreac nspecmx = max(abs(nspec(k)),nspecmx) end do if(nspecmx.gt.maxsp)then stop 'gckwk: nspecmx > maxsp' endif c----------------------------------------------------------------------------- c ... for numerical efficiency do k=1,numreac nk = 0 do ir = 1, nreac(k) nk = nk + 1 kspk(nk,k) = nunk(ir,k) anuk(nk,k) = anu(ir,k) end do do ip3 = 4, nprod(k)+3 nk = nk + 1 kspk(nk,k) = nunk(ip3,k) anuk(nk,k) = anu(ip3,k) end do nkk(k) = nk end do c----------------------------------------------------------------------------- #ifdef VERBOSE open(unit=20,file='stiff.out',status='unknown') write(20,'("# k, nkk(k), kspk(1:6,k), anuk(1:6,k), nkflag")') do k=1,numreac nkflag = 0 do i = 1, nkk(k) kspiv = kspk(i,k) do j = 1, nkk(k) if (i.ne.j .and. kspiv.eq.kspk(j,k))then nkflag = 1 goto 1010 endif end do end do 1010 continue write(20,'(15i4)')k,nkk(k),(kspk(i,k),i=1,6), + (nint(anuk(i,k)),i=1,6),nkflag end do close(20) #endif c----------------------------------------------------------------------------- return end c*****************************************************************************