c ... code to compute jacobian. subroutine dojac(n_slack_species_index,tc,rc,ntmx,temptab, 1 temp1tab,odtab,cgsptab,hsptab,cpsptab, 2 cgsptb_t,hsptb_t,cpsptb_t,cspl, 3 l_rxn,k_spec,p_ropl_p,pd,n_size,pdnc) IMPLICIT NONE include 'parameter.par' include 'chem.com' include 'chem_ref.com' c------------------------------------------------------------------------------- integer ntmx, n_size real*8 temptab(ntmx) real*8 cgsptab(nsmx,ntmx) real*8 hsptab(nsmx,ntmx) real*8 cpsptab(nsmx,ntmx) real*8 cgsptb_t(nsmx,ntmx) real*8 hsptb_t(nsmx,ntmx) real*8 cpsptb_t(nsmx,ntmx) double precision temp1tab,odtab real*8 z(nsmx+2) real*8 zdot(nsmx+2) real*8 cgspl(nsmx),hspl(nsmx),cpspl(nsmx) real*8 yvl(nsmx+3),rr_f(idim),rr_r(idim),aaf(idim),ropl(idim) real*8 qfac(idim),ropl1(idim),cthb(idim),rklow(idim),pr(idim) real*8 prlog(idim),fc(idim),xp(idim),arg(idim),aa(idim) real*8 qff(idim) real*8 bb(idim),clog(idim),oclog(idim),fcent(idim) real*8 yvdot(nsmx+3),yy(nsmx),wspl(nsmx),cspl(nsmx),xik(idim) real*8 econ(5,idim) real*8 p_ropl_p(nsmx+3,idim),pdnc(n_size,n_size),pd(nsmx+1,nsmx+1) real*8 wcheml(nsmx) double precision one,rc integer n_slack_species_index integer i,k,j,kthb,ispthb,ksp,kfl,nk,itest,ntst,neq integer nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v,nfar_v integer npar_v,nrowpd,jtab,jtpl_rxn,k_spec,jtp,l_rxn double precision capsp,gamma,gcon,wbar,tc,otc,vtemp double precision vlntemp,rtab,ctot,vt,vot,hwt,sumw double precision wbarbiro,cp,orcp,ewtl c------------------------------------------------------------------------------- one = 1.0d0 xl10 = log(10.d0) dtc = 0.01 do i=1,ns cspl_diag(i) = ospwt(i) enddo comk oct. '96 do i=1,ns damspwt(i) = spwt(i)*damkohlr enddo c ... identify which reactions have c third-body efficiencies and/or are fall-off reactions c Also, store "constant" part of p_qfac_p into qstore c do k=1,numreac itype(k) = 0 do j=1,ns+3 qstore(j,k) = 0. enddo enddo if(nthb.gt.0)then do kthb=1,nthb k =ithb(kthb) itype(k) = 1 do j=1,ns qstore(j,k) = cspl_diag(j) enddo do ispthb = 1,ntbs(kthb) ksp = nktb(ispthb,kthb) qstore(ksp,k) = qstore(ksp,k) + + (aik(ispthb,kthb) - one)*cspl_diag(ksp) enddo enddo endif if(nfal.gt.0)then do kfl = 1, nfal k = ifal(kfl) itype(k) = 1 enddo endif #ifdef ALL_SPECIES wprims = 0.d0 do i=1,ns wprim(i) = ospwt(i) - wprims enddo #else wprims = ospwt(n_slack_species_index) do i=1,ns wprim(i) = ospwt(i) - wprims enddo wprim(n_slack_species_index) = wprims #endif c c ... generate list of independent species entering c ... in type 0 reactions c do k = 1, numreac if(itype(k).eq.0)then nki(k) = 1 kspi(1,k) = kspk(1,k) do nk = 2,nkk(k) ksp = kspk(nk,k) c c ... test if this is in fact an independent species c itest = 0 do ntst = 1,nki(k) if(ksp.eq.kspi(ntst,k))itest = 1 enddo if(itest.eq.0)then nki(k) = nki(k) + 1 kspi(nki(k),k) = ksp endif enddo endif enddo comk oct. '96 c------------------------------------------------------------------------------- c If we are solving for all the species (ns total) it is c------------------------------------------------------------------------------- c ... z(k) = rho * ys(i,j,k) for k=1,2,...,ns at time t_n c ... z(ns+1) = rho at time t_n c ... z(ns+2) = p0 at time t_n c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c Otherwise if we are solving for ns-1 species it is c------------------------------------------------------------------------------- c ... z(k) = rho * ys(i,j,k) for k=1,2,...,ns-1 at time t_n c ... z(ns) = rho at time t_n c ... z(ns+1) = p0 at time t_n c------------------------------------------------------------------------------- #ifdef ALL_SPECIES do k=1,ns z(k) = cspl(k)*spwt(k) end do z(ns+1) = rc z(ns+2) = p0 neq = ns+2 #else if (n_slack_species_index.eq.1) then do k=1,ns-1 z(k) = cspl(k+1)*spwt(k+1) end do elseif (n_slack_species_index.eq.ns) then do k=1,(ns-1) z(k) = cspl(k)*spwt(k) end do else do k=1,(n_slack_species_index-1) z(k) = cspl(k)*spwt(k) end do do k=n_slack_species_index,ns-1 z(k) = cspl(k+1)*spwt(k+1) end do end if z(ns) = rc z(ns+1) = p0 neq = ns+1 #endif c------------------------------------------------------------------------------- nsmx_v = nsmx idim_v = idim ntmx_v = ntmx maxsp_v = maxsp maxtb_v = maxtb nfar_v = nfar npar_v = npar nrowpd = nsmx+2 c------------------------------------------------------------------------------- capsp = 0.0d0 ! dp0dt c------------------------------------------------------------------------------- gamma = 1.4 gcon = (gamma-1.0)/gamma c------------------------------------------------------------------------------- call gzdot(neq,z,zdot,ns, + nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v,nfar_v,npar_v, + xl10,cspl_diag,dtc, + cgspl,hspl,cpspl,yvl,rr_f,rr_r,aaf,ropl,qfac,ropl1,cthb,rklow,pr, + prlog,fc,xp,arg,aa,qff,bb,clog,oclog,fcent,yvdot,yy,wspl,cspl, + xik,rc,wbar,tc,otc,vtemp,vlntemp,rtab,ctot,vt,vot,capsp,hwt, + sumw,wbarbiro,cp,p0,orcp,econ, + jtab,jtp,tempref,wt,rg,cgref,cpref,href, + ospwt,gcon,rexpo,rcape, + numreac,itfdep,par,cape,nspec,nprod,kflagrev,itrdep, + nkk,anuk,kspk,anusum,nreac,nu,nunk,nthb,ithb, + ntbs,nktb,aik,nfal,ifal,ifop,kfal,capel,fpar,damspwt, + capb,rcapb,capbl,wcheml,damkohlr,ewtl, + temptab,temp1tab,odtab,cgsptab,cpsptab,hsptab, + n_slack_species_index) call jac(neq,pd,n_size,pdnc,nrowpd,ns, + nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v,nfar_v,npar_v, + xl10,cspl_diag,dtc, + cgspl,hspl,cpspl,yvl,rr_f,rr_r,aaf,ropl,qfac,ropl1,cthb,rklow,pr, + prlog,fc,xp,arg,aa,qff,bb,clog,oclog,fcent,yvdot,yy,wspl,cspl, + xik,rc,wbar,tc,otc,vtemp,vlntemp,rtab,ctot,vt,vot,capsp,hwt, + sumw,wbarbiro,cp,p0,orcp,econ, + jtab,jtp,tempref,wt,rg,cgref,cpref,href, + ospwt,gcon,rexpo,rcape, + numreac,itfdep,par,cape,nspec,nprod,kflagrev,itrdep, + nkk,anuk,kspk,anusum,nreac,nu,nunk,nthb,ithb, + ntbs,nktb,aik,nfal,ifal,ifop,kfal,capel,fpar,damspwt, + cgsptb_t,cpsptb_t,hsptb_t, + qstore,wprim,wprims,nki,kspi,itype,p_ropl_p, + nsmx,maxsp,idim,n_slack_species_index,tymref) c------------------------------------------------------------------------------- c write(*,*)' done with dojac' c------------------------------------------------------------------------------- return end