c ... code to compute jacobian. c ... Assumes gzdot has been called before coming here. subroutine jac(neq,pd,n_size,pdnc1,nrowpd,ns, + nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v,nfar_v,npar_v, + xl10,cspl_diag,dtc, + cgspl,hspl,cpspl,yv,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,p0l,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) IMPLICIT NONE c Global variables integer neq,nrowpd,ns,nsmx_v,idim_v,ntmx_v,maxsp_v,maxtb_v integer nfar_v,npar_v,jtab,jtp,numreac,nthb,nfal,nsmx integer maxsp,idim,n_slack_species_index,n_size integer itfdep(idim_v),nspec(idim_v),nprod(idim_v) integer kflagrev(idim_v),itrdep(idim_v),nkk(idim_v) integer kspk(maxsp_v,idim_v),nreac(idim_v) integer nu(maxsp_v,idim_v),nunk(maxsp_v,idim_v) integer ithb(idim_v),ntbs(idim_v),nktb(maxtb_v,idim_v) integer ifal(idim_v),ifop(idim_v),kfal(idim_v) integer nki(idim_v),kspi(maxsp_v,idim_v),itype(idim_v) double precision xl10,dtc,rc,wbar,tc,otc,vtemp,vlntemp,rtab double precision ctot,vt,vot,capsp,hwt,sumw,wbarbiro,cp,p0l double precision orcp,tempref,cgref,cpref,href, tymref double precision gcon,temp1tab,odtab,wprims real*8 pd(nrowpd,nrowpd),pdnc(nsmx+2,nsmx+2),cspl_diag(nsmx_v) real*8 pdnc1(n_size,n_size) real*8 cgspl(nsmx_v),hspl(nsmx_v),cpspl(nsmx_v) real*8 yv(nsmx_v+3),rr_f(idim_v),rr_r(idim_v),aaf(idim_v) real*8 ropl(idim_v),qfac(idim_v),ropl1(idim_v),cthb(idim_v) real*8 rklow(idim_v),pr(idim_v),prlog(idim_v),fc(idim_v) real*8 xp(idim_v),arg(idim_v),aa(idim_v),qff(idim_v) real*8 bb(idim_v),clog(idim_v),oclog(idim_v),fcent(idim_v) real*8 yvdot(nsmx_v+3),yy(nsmx_v),wspl(nsmx_v),cspl(nsmx_v) real*8 xik(idim_v),econ(5,idim_v),wt(nsmx_v),rg(nsmx_v) real*8 ospwt(nsmx),rexpo(idim_v),rcape(idim_v) real*8 par(npar_v,idim_v),cape(idim_v),anuk(maxsp_v,idim_v) real*8 anusum(idim_v),aik(maxtb_v,idim_v),capel(idim_v) real*8 fpar(nfar_v,idim_v),damspwt(nsmx_v) real*8 cgsptb_t(nsmx_v,ntmx_v),cpsptb_t(nsmx_v,ntmx_v) real*8 hsptb_t(nsmx_v,ntmx_v),temptab(ntmx_v) real*8 qstore(nsmx_v+3,idim_v),wprim(nsmx_v) real*8 p_ropl_p(nsmx+3,idim) c Local Variables integer k,j,nk integer isp,nyp,ispk,nid,ip3,ip3k,kthb integer kfl,kfc,ifp,i,ksp,nk2 double precision cgsv,cpv,hsav,rrfkl,rrrkl double precision rrr_t,xik_t,crprod,hsv,cvp,crp2,crp1 double precision cpprod,cpp2,cpp1,rklow_t,p_xp_p_pr double precision comm,fcent_t,gg,af,ap,bf,bp double precision gcom,p_fc_p_fcent,p_fc_p_pr,pr1,prcm,der double precision pcor,ankl,wb2,rocp double precision orcp2,cpw,sum,sum1,sum2 double precision sum3,sum4,sum5,sum6 double precision rowbar,rowb2,p_t_p_rho,p_t_p_p0,rcsq double precision sum_a,sum_b,sum_c,sum_d,dyndyj,dyndt,wsum double precision pots2,drdt,drdyj real*8 cgspl_T(nsmx),cgdif(nsmx) real*8 cpspl_T(nsmx),cpdif(nsmx) real*8 hspl_T(nsmx),hsdif(nsmx) real*8 rrf_T(idim) real*8 crpvec(maxsp),cppvec(maxsp) real*8 crprod_p(nsmx+3),cpprod_p(nsmx+3) real*8 p_ropl_old_p(nsmx+3,idim) real*8 p_qfac_p(nsmx+3,idim) real*8 p_cthb_p(nsmx+3) real*8 p_pr_p(nsmx+3) real*8 p_fc_p(nsmx+3) real*8 p_pcor_p(nsmx+3) real*8 p_wb_p(nsmx+3) real*8 p_cpw_p(nsmx+3) real*8 p_yns_p(nsmx+3) real*8 p_T_p(nsmx+3) real*8 sum_ypdf(nsmx) real*8 pdf(nsmx+3,nsmx+3),pdfnc(nsmx+3,nsmx+3) c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... get derivative w.r.t. T of gibbs free energies c ... get derivative w.r.t. T of c_{p,i} c ... get derivative w.r.t. T of enthalpies c------------------------------------------------------------------------------- c ... linear interpolation from derivative tables do k=1,ns cgsv = cgsptb_t(k,jtab) cgspl_T(k) = cgsv + rtab*(cgsptb_t(k,jtp) - cgsv) cpv = cpsptb_t(k,jtab) cpspl_T(k) = cpv + rtab*(cpsptb_t(k,jtp) - cpv) hsv = hsptb_t(k,jtab) hspl_T(k) = hsv + rtab*(hsptb_t(k,jtp) - hsv) end do c------------------------------------------------------------------------------ c c c species Source Terms c c c------------------------------------------------------------------------------ do k = 1, numreac comk oct. '96 c ... this loop has been moved from below do j=1,ns+1 p_ropl_p(j,k) = 0. end do c ... statements have been moved from below rrfkl = rr_f(k) rrrkl = rr_r(k) comk oct. '96 c ... evaluate forward reaction rate, rr_f if(itfdep(k).ne.0)then rrf_T(k) = rr_f(k)*(par(2,k)+cape(k)*otc)*otc else rrf_T(k) = 0. endif c------------------------------------------------------------------------------- c ... nspec(k) > 0 iff reaction is reversible c ... xik = SUM_i=1,ns nu(i,k)*G_ik c ... evaluate (1/equilibrium_reaction_rate) = orr_e c ... evaluate reverse reaction rate, rr_r c ... kflagrev(k) = 1 -> reverse Arrhenius rate parameters are defined c ... use these rather than orr_e to find rr_r if(nspec(k).gt.0)then if(kflagrev(k).eq.1)then if(itrdep(k).ne.0)then rrr_T = rr_r(k)*(rexpo(k)+rcape(k)*otc)*otc else rrr_T = 0. endif else xik_T = 0.0 do nk = 1, nkk(k) xik_T = xik_T + anuk(nk,k)*cgspl_T(kspk(nk,k)) end do rrr_T = rrf_T(k)*aaf(k) + + rr_r(k)*(xik_T - xik(k)*otc + anusum(k))*otc endif else rrr_T = 0.0 endif c------------------------------------------------------------------------------- c ... find reactants concentrations product PROD(c_j**nu_jk^prime)_j=1,ns c ... it is assumed below that nyp is a positive integer (non-zero) c ... which should be true in general crprod = 1.0 do j=1,maxsp crpvec(j) = 1.0 end do comk oct. '96 c ... I don't need this loop anymore c do j=1,ns+3 c crprod_p(j) = 0. c end do comk oct. '96 do isp = 1,nreac(k) nyp = - nu(isp,k) cvp = cspl(nunk(isp,k)) crp2 = cvp**(nyp-1) crp1 = crp2 * cvp crprod = crprod*crp1 do ispk=1,nreac(k) if(ispk.eq.isp)then crpvec(ispk) = dble(nyp) * crpvec(ispk) * crp2 else crpvec(ispk) = crpvec(ispk) * crp1 endif end do end do do isp = 1,nreac(k) nid = nunk(isp,k) comk oct. '96 c ... statement eliminated and result directly substituted c into p_ropl_p array c crprod_p(nid) = crpvec(isp) * cspl_diag(nid) c p_ropl_p(nid,k) = rrfkl * crpvec(isp) * cspl_diag(nid) comk oct. '96 end do c------------------------------------------------------------------------------- c ... for reverse reaction, c ... find products concentrations product PROD(c_j**nu_jk^double-prime)_j=1,ns c ... nspec(k) > 0 if reaction is reversible c ... it is assumed below that nu(ip3,k) is a positive integer (non-zero) c ... which should be true in general if(nspec(k).gt.0)then cpprod = 1.0 do j=1,maxsp cppvec(j) = 1.0 end do comk oct. '96 c ... I don't need this loop anymore c do j=1,ns+3 c cpprod_p(j) = 0. c end do comk oct. '96 do ip3 = 4,nprod(k)+3 nyp = nu(ip3,k) cvp = cspl(nunk(ip3,k)) cpp2 = cvp**(nyp-1) cpp1 = cpp2 * cvp cpprod = cpprod * cpp1 do ip3k = 4,nprod(k)+3 if(ip3k.eq.ip3)then cppvec(ip3k) = dble(nyp) * cppvec(ip3k) * cpp2 else cppvec(ip3k) = cppvec(ip3k) * cpp1 endif end do end do do ip3 = 4,nprod(k)+3 nid = nunk(ip3,k) comk oct. '96 c ... statement eliminated and result directly substituted c into p_ropl_p array c cpprod_p(nid) = cppvec(ip3) * cspl_diag(nid) c p_ropl_p(nid,k) = p_ropl_p(nid,k) - 1 rrrkl * cppvec(ip3) * cspl_diag(nid) comk oct. '96 end do else cpprod = 0. comk oct. '96 c ... I don't need this loop anymore c do j=1,ns+3 c cpprod_p(j) = 0. c end do endif c------------------------------------------------------------------------------- c ... evaluate reaction rate-of-progress comk oct. '96 c ... this loop has been moved up c do j=1,ns+3 c p_ropl_p(j,k) = 0. c end do comk oct. '96 p_ropl_p(ns+1,k) = rrf_T(k)*crprod - rrr_T*cpprod comk oct. '96 c ... these statements have been moved to beginning of loop c rrfkl = rr_f(k) c rrrkl = rr_r(k) comk oct. '96 comk oct. '96 c ... this loop has been eliminated! c ... the operations have been moved into c the calculation of product terms above c in order to avoid filling (needlessly) c zero entries c c do j=1,ns c p_ropl_p(j,k) = rrfkl * crprod_p(j) - rrrkl * cpprod_p(j) c end do comk oct. '96 c------------------------------------------------------------------------------- end do c ... end of loop over all reactions for rate-of-progress c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... initialize qfac array, product of concentrations to be used in c ... third body and fall off reactions do k=1,numreac comk oct. '96 c ... save some operations if the reaction c has neither third body efficiencies nor c is fall-off if(itype(k).ne.0)then do j=1,ns+1 p_qfac_p(j,k) = 0. end do endif end do c ... evaluate third-body-scaled reaction rates-of-progress c ... nb. qfac is the product of concentrations, with possible 3rd c ... body efficiencies. c ... if the reaction is also a fall-off reaction with kfal=0, then c ... qfac is just used to find pr below, and then reset to 1.0 c ... otherwise, it is actually used to scale the rate of progress c ... later on below if(nthb.gt.0)then do kthb = 1,nthb k = ithb(kthb) comk oct. '96 c ... this loop has been moved to otherinit c do j=1,ns c p_qfac_p(j,k) = cspl_diag(j) c end do comk oct. '96 comk oct. '96 c ... this loop has been moved to otherinit c do ispthb = 1,ntbs(kthb) c ksp = nktb(ispthb,kthb) c p_qfac_p(ksp,k) = p_qfac_p(ksp,k) + c + (aik(ispthb,kthb)-1.0d0)*cspl_diag(ksp) c end do comk oct. '96 c ... load p_qfac_p from qstore c which is computed in otherinit and saved in common prepjac do j=1,ns p_qfac_p(j,k) = qstore(j,k) enddo comk oct. '96 comk oct. '96 c ... the following two loops have been blended into one c c do j=1,ns+3 c p_ropl_old_p(j,k) = p_ropl_p(j,k) c end do c c do j=1,ns+3 c p_ropl_p(j,k) = p_qfac_p(j,k)*ropl1(k) + c + qff(k)*p_ropl_p(j,k) c end do do j=1,ns+1 p_ropl_old_p(j,k) = p_ropl_p(j,k) p_ropl_p(j,k) = p_qfac_p(j,k)*ropl1(k) + + qff(k)*p_ropl_p(j,k) end do comk oct. '96 end do endif c------------------------------------------------------------------------------- c ... evaluate fall-off pressure corrections c ... vt is dimensional temperature (K), vot = 1/vt c ... these are used in the computation of the dimensionless F if(nfal.gt.0)then do kfl = 1, nfal k = ifal(kfl) kfc = kfal(kfl) ifp = ifop(kfl) if(kfc.eq.0)then do j=1,ns+1 p_ropl_p(j,k) = p_ropl_old_p(j,k) p_cthb_p(j) = p_qfac_p(j,k) p_qfac_p(j,k) = 0. end do else do j=1,ns+1 p_cthb_p(j) = 0. end do p_cthb_p(kfc) = cspl_diag(kfc) endif rklow_T = rklow(k)*(fpar(2,kfl)+capel(k)*otc)*otc do j=1,ns+1 p_pr_p(j) = rklow(k) * p_cthb_p(j) / rr_f(k) end do j = ns+1 p_pr_p(j) = p_pr_p(j) + rklow_T * cthb(k) / rr_f(k) + - rklow(k) * cthb(k) * rrf_T(k) / (rr_f(k)*rr_f(k)) if(ifp.eq.1)then ! Lindemann form do j=1,ns+1 p_fc_p(j) = 0. end do else ! ifop = 2, 3, or 4 if(ifp.eq.2)then ! SRI form p_xp_p_pr = -2.*prlog(k)*xp(k)*xp(k)/pr(k)/xl10 ! O28 do j=1,ns+1 p_fc_p(j) = 0. end do comm = fc(k)*log(arg(k))*p_xp_p_pr do j=1,ns+1 p_fc_p(j) = comm*p_pr_p(j) ! O29 end do j=ns+1 p_fc_p(j) = p_fc_p(j) + fc(k) * ( + otc*fpar(8,kfl) + xp(k) * ( + econ(1,k)*fpar(5,kfl)*vot*otc - + (tempref/fpar(6,kfl))*econ(2,k) + )/arg(k) ) ! O27 else ! 6 or 7-param TROE form fcent_T = - econ(3,k)/fpar(5,kfl) - + econ(4,k)/fpar(6,kfl) fcent_T = fcent_T*tempref ! O13 if(ifp.eq.4)then ! 7-parameter TROE form fcent_T = fcent_T + fpar(7,kfl)*vot*otc* + econ(5,k) ! O13 endif gg = 2.*oclog(k)*oclog(k)*clog(k)/(bb(k)*bb(k)) ! O15 af = -0.67/(xl10*fcent(k)) ! O17 ap = 1./(xl10*pr(k)) ! O18 bf = -1.1762/(xl10*fcent(k)) ! O19 bp = -0.14/(xl10*pr(k)) ! O20 gcom = fc(k)*log(fcent(k))*gg p_fc_p_fcent = oclog(k)*fc(k)/fcent(k) - + gcom*(af*bb(k) - aa(k)*bf) ! O14 p_fc_p_pr = - gcom*(ap*bb(k) - aa(k)*bp) ! O16 comk oct. '96 c ... unroll the loop c do j=1,ns+3 c p_fc_p(j) = 0. c end do comk oct. '96 do j=1,ns+1 p_fc_p(j) = p_fc_p_pr * p_pr_p(j) ! O22 end do j=ns+1 p_fc_p(j) = p_fc_p(j) + p_fc_p_fcent * fcent_T ! O21 comk oct. '96 c p_fc_p(ns+2) = 0. c p_fc_p(ns+3) = 0. comk oct. '96 endif endif pr1 = 1.0+pr(k) prcm = pr(k)/pr1 der = fc(k)/(pr1*pr1) pcor = fc(k)*prcm do j=1,ns+1 p_pcor_p(j) = p_fc_p(j) * prcm + der*p_pr_p(j) end do do j=1,ns+1 p_ropl_p(j,k) = p_ropl_p(j,k)*pcor + + ropl(k)*p_pcor_p(j)/pcor end do end do endif c------------------------------------------------------------------------------- c ... evaluate contributions to reactants and products reaction rates c c ... initialize pdf c do j=1,ns+3 do i=1,ns pdf(i,j) = 0. end do end do do k=1,numreac comk oct. '96 c ... load only non-zero entries c if(itype(k).eq.0)then c ... for type 0 reactions, p_ropl_p has non-zero c entries only when the derivatives is w.r.t. c the corresponding species concentrations or c temperature do nk = 1, nkk(k) ksp = kspk(nk,k) ankl = anuk(nk,k) do nk2=1,nki(k) j = kspi(nk2,k) pdf(ksp,j) = pdf(ksp,j) + ankl*p_ropl_p(j,k) end do pdf(ksp,ns+1) = pdf(ksp,ns+1) + ankl*p_ropl_p(ns+1,k) end do else do nk = 1, nkk(k) ksp = kspk(nk,k) ankl = anuk(nk,k) do j=1,ns+1 pdf(ksp,j) = pdf(ksp,j) + ankl*p_ropl_p(j,k) end do end do endif comk oct. '96 end do c------------------------------------------------------------------------------- c ... final scaling comk oct. '96 c ... take advantage of damspwt in order to save (ns)^2 multiplications c do j=1,ns+3 c do ksp=1,ns c pdf(ksp,j) = pdf(ksp,j)*spwt(ksp)*damkohlr c end do c end do do j=1,ns+1 do ksp=1,ns pdf(ksp,j) = pdf(ksp,j)*damspwt(ksp) end do end do comk oct. '96 c c ................................................................. c c Density Source Term c c ................................................................. c c ... derivative of the source term w.r.t. c ... expanded state vector! c ... first compute some useful expressions wb2 = wbarbiro*wbarbiro do j=1,ns p_wb_p(j) = -wb2*ospwt(j) ! O47 end do rocp = rc*cp orcp2 = orcp*orcp cpw = cp*wbar do j=1,ns p_cpw_p(j) = wbarbiro*cpspl(j) + rocp*p_wb_p(j) ! O49 end do sum = 0. do i=1,ns sum = sum + yv(i)*cpspl_T(i) end do p_cpw_p(ns+1) = wbarbiro*sum ! O50 c p_cpw_p(ns+2) = 0. ! O51 c p_cpw_p(ns+3) = 0. ! O52 do j=1,ns+3 pdf(ns+2,j) = 0. end do do j=1,ns sum1 = 0. sum2 = 0. do i=1,ns sum1 = sum1 + hspl(i)*pdf(i,j) sum2 = sum2 + pdf(i,j)*ospwt(i) end do pdf(ns+2,j) = orcp*rc/tc*sum1 + - cpspl(j)*orcp2*rc/tc*hwt + + rc*capsp/p0l*gcon/cpw/cpw*p_cpw_p(j) + - rc*p_wb_p(j)*sumw + - rc*wbarbiro*sum2 ! O41 end do j = ns+1 sum3 = 0. sum4 = 0. sum5 = 0. sum6 = 0. do i=1,ns sum3 = sum3 + yvdot(i)*hspl_T(i) sum4 = sum4 + hspl(i)*pdf(i,j) sum5 = sum5 + pdf(i,j)*ospwt(i) sum6 = sum6 + yv(i)*cpspl_T(i) end do pdf(ns+2,j) = -orcp*rc/tc/tc*hwt + + orcp*rc/tc*sum3 + + orcp*rc/tc*sum4 + - orcp2*rc/tc*sum6*hwt + + rc*capsp/p0l*gcon/cpw/cpw*p_cpw_p(j) + - rc*wbarbiro*sum5 ! O42 j = ns+2 pdf(ns+2,j) = orcp*hwt/tc + + capsp/p0l*(1. - gcon/cpw) + - wbarbiro*sumw ! O43 j = ns+3 pdf(ns+2,j) = -rc*capsp/p0l/p0l*(1. - gcon/cpw) ! O44 c c ................................................................. c c Pressure c c ................................................................. c c c do j=1,ns+3 c pdf(ns+3,j) = 0. c end do c c ................................................................. c c Reduced State Vector c c ................................................................. c c ... zdot first c c ... concentration of species "ns" and temperature c ... have been eliminated from the state vector!! c c ... The ``trick'' now is to reduce the Jacobian pdf --> pd c ... Chain Rule Galore... comk oct. '96 c ... We don't actually need these c do j=1,ns-1 c p_yns_p(j) = -1. c end do c p_yns_p_rho = 1. c wprims = ospwt(ns) c do i=1,ns-1 c wprim(i) = ospwt(i) - wprims c end do rowbar = rc/wbar rowb2 = rowbar*rowbar #ifdef ALL_SPECIES do j=1,ns p_T_p(j) = - p0l*wprim(j)/rowb2 ! O38 end do #else if (n_slack_species_index.eq.1) then do j=1,ns-1 p_T_p(j) = - p0l*wprim(j+1)/rowb2 ! O38 end do elseif (n_slack_species_index.eq.ns) then do j=1,ns-1 p_T_p(j) = - p0l*wprim(j)/rowb2 ! O38 end do else do j=1,(n_slack_species_index-1) p_T_p(j) = - p0l*wprim(j)/rowb2 enddo do j=n_slack_species_index,ns-1 p_T_p(j) = - p0l*wprim(j+1)/rowb2 enddo endif #endif p_T_p_rho = - p0l*wprims/rowb2 ! O37 p_T_p_p0 = 1./rowbar ! O36 comk oct. '96 c ... substitute the value of p_yns_p(j) = -1. c #ifdef ALL_SPECIES do j=1,ns do i=1,ns pd(i,j) = pdf(i,j) + + pdf(i,ns+1)*p_T_p(j) end do end do #else if (n_slack_species_index.eq.1) then do j=2,ns do i=2,ns pd(i-1,j-1) = pdf(i,j) - + pdf(i,n_slack_species_index) + + pdf(i,ns+1)*p_T_p(j) enddo enddo elseif (n_slack_species_index.eq.ns) then do j=1,ns-1 do i=1,ns-1 pd(i,j) = pdf(i,j) - + pdf(i,n_slack_species_index) + + pdf(i,ns+1)*p_T_p(j) end do end do else do j=1,(n_slack_species_index-1) do i=1,(n_slack_species_index-1) pd(i,j) = pdf(i,j) - + pdf(i,n_slack_species_index) + + pdf(i,ns+1)*p_T_p(j) enddo enddo do j=n_slack_species_index+1,ns do i=n_slack_species_index+1,ns pd(i-1,j-1) = pdf(i,j) - + pdf(i,n_slack_species_index) + + pdf(i,ns+1)*p_T_p(j) enddo enddo endif #endif comk oct. '96 comk oct. '96 c ... substitute the values of p_yns_p_rho = 1. c #ifdef ALL_SPECIES do i=1,ns pd(i,ns+1) = pdf(i,ns+1)*p_T_p_rho pd(i,ns+2) = pdf(i,ns+1)*p_T_p_p0 end do comk oct. '96 comk oct. '96 c ... substitute the value of p_yns_p(j) = -1. c do j=1,ns pd(ns+1,j) = pdf(ns+2,j) + + pdf(ns+2,ns+1)*p_T_p(j) end do comk oct. '96 comk oct. '96 c ... substitute the values of p_yns_p_rho = 1. c pd(ns+1,ns+1) = pdf(ns+2,ns+2) + + pdf(ns+2,ns+1)*p_T_p_rho comk oct. '96 pd(ns+1,ns+2) = pdf(ns+2,ns+3) + pdf(ns+2,ns+1)*p_T_p_p0 do j=1,ns+2 pd(ns+2,j) = 0. end do #else do i=1,ns-1 c pd(i,ns) = pdf(i,ns+2) + pdf(i,ns)*p_yns_p_rho + c + pdf(i,ns+1)*p_T_p_rho c pd(i,ns+1) = pdf(i,ns+3) + pdf(i,ns+1)*p_T_p_p0 pd(i,ns) = pdf(i,n_slack_species_index) + + pdf(i,ns+1)*p_T_p_rho pd(i,ns+1) = pdf(i,ns+1)*p_T_p_p0 end do comk oct. '96 comk oct. '96 c ... substitute the value of p_yns_p(j) = -1. c do j=1,ns-1 c pd(ns,j) = pdf(ns+2,j) + pdf(ns+2,ns)*p_yns_p(j) + c + pdf(ns+2,ns+1)*p_T_p(j) pd(ns,j) = pdf(ns+2,j) - pdf(ns+2,n_slack_species_index) + + pdf(ns+2,ns+1)*p_T_p(j) end do comk oct. '96 comk oct. '96 c ... substitute the values of p_yns_p_rho = 1. c c pd(ns,ns) = pdf(ns+2,ns+2) + pdf(ns+2,ns)*p_yns_p_rho + c + pdf(ns+2,ns+1)*p_T_p_rho pd(ns,ns) = pdf(ns+2,ns+2) + pdf(ns+2,n_slack_species_index) + + pdf(ns+2,ns+1)*p_T_p_rho comk oct. '96 pd(ns,ns+1) = pdf(ns+2,ns+3) + pdf(ns+2,ns+1)*p_T_p_p0 do j=1,ns+1 pd(ns+1,j) = 0. end do #endif c-------------------------------------------------------------------------------- c ... with N = number of species = ns, c ... the jacobian pdf (N+3)x(N+3) is the jacobian of the following full ODE system: c ... c ... du/dt = f(u) where u is (N+3)x1, u = {rho*Y_1,...,rho*Y_N,T,rho,P} c ... and f(u)={Da*w_1,...,Da*w_N,S_T,S_rho,S_P} c ... where source terms S_T, S_rho are too messy to put here,j and S_P=0 c ... c ... if we are solving for (ns-1) species it is: c ... pd (N+1)x(N+1) is the jacobian of the reduced ODE system : c ... c ... dv/dt = g(v) where v is (N+1)x1, v = {rho*Y_1,...,rho*Y_[N-1],rho,P} c ... and g(v)={Da*w_1,...,Da*w_[N-1],S_rho,S_P} c ... (the reduced system makes use of the state equation & Sum(Y_i)=1.0) c ... otherwise if we are solving for all the species (ns) it is: c ... pd (N+2)x(N+2) is the jacobian of the reduced ODE system : c ... c ... dv/dt = g(v) where v is (N+2)x1, v = {rho*Y_1,...,rho*Y_[N],rho,P} c ... and g(v)={Da*w_1,...,Da*w_[N],S_rho,S_P} c ... (the reduced system makes use of the state equation) c-------------------------------------------------------------------------------- c ... we now fill out the pdfnc and pdnc jacobians for non-dflame uses c ... the nc ending implies non-conservative c ... c ... the jacobian pdfnc (N+3)x(N+3) is the jacobian of the full ODE system: c ... c ... du/dt = f(u) where u is (N+3)x1, u = {Y_1,...,Y_N,T,rho,P} c ... and f(u)={Da*w_1/rho,...,Da*w_N/rho,S_T,S_rho,S_P} c ... c ... if we are solving for (ns-1) species it is: c ... then pdnc (N+1)x(N+1) is the jacobian of the reduced ODE system: c ... c ... dv/dt = g(v) where v is (N+1)x1, v = {Y_1,...,Y_[N-1],T,P} c ... and g(v)={Da*w_1/rho,...,Da*w_[N-1]/rho,S_T,S_P} c ... c ... otherwise if we are solving for all the species (ns) it is: c ... then pdnc (N+2)x(N+2) is the jacobian of the reduced ODE system: c ... c ... dv/dt = g(v) where v is (N+2)x1, v = {Y_1,...,Y_[N],T,P} c ... and g(v)={Da*w_1/rho,...,Da*w_[N]/rho,S_T,S_P} c ... So, we proceed to first fill up pdfnc, and then pdnc c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- c ... pdfnc c ... c ... fill up pdfnc(i,j) with i=1,ns+1 j=1,ns+2 ! tc = T ! rc = rho rcsq = rc*rc ! i=1,ns, j=1,ns subset of pdfnc do j=1,ns do i=1,ns pdfnc(i,j) = pdf(i,j) end do end do ! i=1,ns j=ns+1 do i=1,ns pdfnc(i,ns+1) = pdf(i,ns+1)/rc end do ! i=1,ns, j=ns+2 ! sum_ypdf(i) = sum_{k=1}^ns (rho*Y_k) * pdf(i,k) ! yvdot(i) = Da w_i, i=1,2,...,ns do i=1,ns sum_ypdf(i) = 0.0d0 do k=1,ns sum_ypdf(i) = sum_ypdf(i) + yv(k) * pdf(i,k) end do end do do i=1,ns pdfnc(i,ns+2) = (sum_ypdf(i) - yvdot(i) ) / rcsq end do ! i=ns+1 j=1,ns ! hwt = sum_k=1,..,ns { Da h_k w_k } ! cpspl(j) = cp_j ! cp = cp_mixture ! orcp = 1/(rho*cp) do j=1,ns sum1 = 0.0d0 do k=1,ns sum1 = sum1 + hspl(k) * pdf(k,j) end do pdfnc(ns+1,j) = (-sum1 + cpspl(j)*orcp*hwt)/cp end do ! i=ns+1 j=ns+1 ! sum_a = sum_k=1,ns {Da w_k dh_k/dT} ! sum_b = sum_k=1,ns {h_k pdf(k,ns+1)} ! sum_c = {1/rho}*sum_k=1,ns {rho*Y_k dcp_k/dT} = d cp / d T ! = sum_k=1,ns {Y_k dcp_k/dT} = d cp / d T sum_a = 0.0d0 sum_b = 0.0d0 sum_c = 0.0d0 do i=1,ns sum_a = sum_a + yvdot(i)*hspl_T(i) sum_b = sum_b + hspl(i)*pdf(i,ns+1) sum_c = sum_c + yv(i)*cpspl_T(i) end do sum_c = sum_c / rc pdfnc(ns+1,ns+1) = orcp * (- sum_a - sum_b + sum_c*hwt/cp) ! i=ns+1 j=ns+2 sum_d = 0.0d0 do k=1,ns sum_d = sum_d + hspl(k)*sum_ypdf(k) end do pdfnc(ns+1,ns+2) = orcp * (-sum_d + hwt) / rc c-------------------------------------------------------------------------------- c-------------------------------------------------------------------------------- c ... pdnc c ... fill up pdnc(i,j) i=1,ns j=1,ns dyndyj = -1.0d0 dyndt = 0.0d0 sum = 0.0d0 #ifdef ALL_SPECIES do i=1,ns sum = sum + (yv(i)/rc) * wprim(i) end do #else if (n_slack_species_index.eq.1) then do i=1,ns-1 sum = sum + (yv(i+1)/rc) * wprim(i+1) end do elseif (n_slack_species_index.eq.ns) then do i=1,ns-1 sum = sum + (yv(i)/rc) * wprim(i) end do else do i=1,(n_slack_species_index-1) sum = sum + (yv(i)/rc) * wprim(i) enddo do i=n_slack_species_index,ns-1 sum = sum + (yv(i+1)/rc) * wprim(i+1) enddo endif #endif wsum = sum + wprims pots2 = -(p0l/tc) / (wsum*wsum) drdt = - (p0l/(tc*tc))/wsum c=============================================================== #ifdef ALL_SPECIES ! i=1,ns j=1,ns do j=1,ns drdyj = pots2 * wprim(j) do i=1,ns pdnc(i,j)= pdfnc(i,j) & + pdfnc(i,ns+2) * drdyj end do end do #else ! i=1,ns-1 j=1,ns-1 if (n_slack_species_index.eq.1) then do j=1,ns-1 drdyj = pots2 * wprim(j+1) do i=1,ns-1 pdnc(i,j)= pdfnc(i+1,j+1) & + pdfnc(i+1,n_slack_species_index) * dyndyj & + pdfnc(i+1,ns+2) * drdyj end do end do elseif (n_slack_species_index.eq.ns) then do j=1,ns-1 drdyj = pots2 * wprim(j) do i=1,ns-1 pdnc(i,j)= pdfnc(i,j) & + pdfnc(i,ns) * dyndyj & + pdfnc(i,ns+2) * drdyj end do end do else do j=1,n_slack_species_index-1 drdyj = pots2 * wprim(j) do i=1,n_slack_species_index-1 pdnc(i,j)= pdfnc(i,j) & + pdfnc(i,n_slack_species_index) * dyndyj & + pdfnc(i,ns+2) * drdyj end do do i=n_slack_species_index,ns-1 pdnc(i,j)= pdfnc(i+1,j) & + pdfnc(i+1,n_slack_species_index) * dyndyj & + pdfnc(i+1,ns+2) * drdyj end do end do do j=n_slack_species_index,ns-1 drdyj = pots2 * wprim(j+1) do i=n_slack_species_index,ns-1 pdnc(i,j)= pdfnc(i+1,j+1) & + pdfnc(i+1,n_slack_species_index) * dyndyj & + pdfnc(i+1,ns+2) * drdyj end do do i=1,n_slack_species_index-1 pdnc(i,j)= pdfnc(i,j+1) & + pdfnc(i,n_slack_species_index) * dyndyj & + pdfnc(i,ns+2) * drdyj end do end do endif #endif c============================================================ #ifdef ALL_SPECIES ! i=1,ns j=ns+1 do i=1,ns pdnc(i,ns+1) = pdfnc(i,ns+1) & + pdfnc(i,ns+2) * drdt end do #else ! i=1,ns-1 j=ns if (n_slack_species_index.eq.1) then do i=1,ns-1 pdnc(i,ns) = pdfnc(i+1,ns+1) & + pdfnc(i+1,n_slack_species_index) * dyndt & + pdfnc(i+1,ns+2) * drdt end do elseif (n_slack_species_index.eq.ns) then do i=1,ns-1 pdnc(i,ns) = pdfnc(i,ns+1) & + pdfnc(i,n_slack_species_index) * dyndt & + pdfnc(i,ns+2) * drdt end do else do j=1,n_slack_species_index-1 pdnc(i,ns) = pdfnc(i,ns+1) & + pdfnc(i,n_slack_species_index) * dyndt & + pdfnc(i,ns+2) * drdt end do do j=n_slack_species_index,ns-1 pdnc(i,ns) = pdfnc(i+1,ns+1) & + pdfnc(i+1,n_slack_species_index) * dyndt & + pdfnc(i+1,ns+2) * drdt enddo endif #endif c============================================================ #ifdef ALL_SPECIES ! i=ns+1 j=1,ns do j=1,ns drdyj = pots2 * wprim(j) pdnc(ns+1,j) = pdfnc(ns+1,j) & + pdfnc(ns+1,ns+2) * drdyj end do #else ! i=ns j=1,ns-1 if (n_slack_species_index.eq.1) then do j=1,ns-1 drdyj = pots2 * wprim(j+1) pdnc(ns,j) = pdfnc(ns+1,j+1) & + pdfnc(ns+1,n_slack_species_index) * dyndyj & + pdfnc(ns+1,ns+2) * drdyj end do elseif (n_slack_species_index.eq.ns) then do j=1,ns-1 drdyj = pots2 * wprim(j) pdnc(ns,j) = pdfnc(ns+1,j) & + pdfnc(ns+1,n_slack_species_index) * dyndyj & + pdfnc(ns+1,ns+2) * drdyj end do else do j=1,n_slack_species_index-1 drdyj = pots2 * wprim(j) pdnc(ns,j) = pdfnc(ns+1,j) & + pdfnc(ns+1,n_slack_species_index) * dyndyj & + pdfnc(ns+1,ns+2) * drdyj end do do j=n_slack_species_index,ns-1 drdyj = pots2 * wprim(j+1) pdnc(ns,j) = pdfnc(ns+1,j+1) & + pdfnc(ns+1,n_slack_species_index) * dyndyj & + pdfnc(ns+1,ns+2) * drdyj end do endif #endif c============================================================ #ifdef ALL_SPECIES c i=ns+1, j=ns+1 pdnc(ns+1,ns+1) = pdfnc(ns+1,ns+1) & + pdfnc(ns+1,ns+2) * drdt #else c i=ns j=ns pdnc(ns,ns) = pdfnc(ns+1,ns+1) & + pdfnc(ns+1,n_slack_species_index) * dyndt & + pdfnc(ns+1,ns+2) * drdt #endif c============================================================ c-------------------------------------------------------------------------------- c Final scaling to dimensionalize #ifdef ALL_SPECIES pdnc1(1,1) = pdnc(ns+1,ns+1) / tymref #ifdef NON_DIM do j=2,(ns+1) pdnc1(1,j) = pdnc(ns+1,j-1 ) enddo do i=2,(ns+1) pdnc1(i,1) = pdnc(i-1,ns+1) enddo do i=2,(ns+1) do j=2,(ns+1) pdnc1(i,j) = pdnc(i-1,j-1) enddo enddo #else do j=2,(ns+1) pdnc1(1,j) = pdnc(ns+1,j-1) * tempref / tymref enddo do i=2,(ns+1) pdnc1(i,1) = pdnc(i-1,ns+1) / (tempref * tymref) enddo do i=2,(ns+1) do j=2,(ns+1) pdnc1(i,j) = pdnc(i-1,j-1) / tymref enddo enddo #endif #else pdnc1(1,1) = pdnc(ns,ns) / tymref #ifdef NON_DIM do j=2,ns pdnc1(1,j) = pdnc(ns,j-1) enddo do i=2,ns pdnc1(i,1) = pdnc(i-1,ns) enddo do i=2,ns do j=2,ns pdnc1(i,j) = pdnc(i-1,j-1) enddo enddo #else do j=2,ns pdnc1(1,j) = pdnc(ns,j-1) * tempref / tymref enddo do i=2,ns pdnc1(i,1) = pdnc(i-1,ns) / (tempref * tymref) enddo do i=2,ns do j=2,ns pdnc1(i,j) = pdnc(i-1,j-1) / tymref enddo enddo #endif #endif return end