c c ========================================================== subroutine flg1eurhok(q,mx,lb,ub,qf,mxo,lbo,ubo, & lbr,ubr,shaper,meqn,mcomp,time) c ========================================================== c c # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c c # Copyright (C) 2003-2007 California Institute of Technology c # Ralf Deiterding, ralf@cacr.caltech.edu c implicit double precision(a-h,o-z) include "ck.i" c integer meqn, mx, mxo dimension q(meqn,mx), qf(mxo) c integer lb(1), ub(1), lbo(1), ubo(1), lbr(1), ubr(1), shaper(1), & mresult, stride, imin(1), imax(1), i, getindx c stride = (ub(1) - lb(1))/(mx-1) imin(1) = max(lb(1), lbo(1), lbr(1)) imax(1) = min(ub(1), ubo(1), ubr(1)) if (mod(imin(1)-lb(1),stride) .ne. 0) then imin(1) = imin(1) + stride - mod(imin(1)-lb(1),stride) endif imin(1) = getindx(imin(1), lb(1), stride) if (mod(imax(1)-lb(1),stride) .ne. 0) then imax(1) = imax(1) - mod(imax(1)-lb(1),stride) endif imax(1) = getindx(imax(1), lb(1), stride) do 10 i = imin(1), imax(1) c Compute rho if (mcomp.eq.1) then qf(i) = 0.d0 do k = 1, Nsp qf(i) = qf(i) + q(k,i) enddo c # Convert g/cm**3 into kg/m**3 if (ckunits) qf(i) = qf(i)*1.d3 c Compute p else if (mcomp.eq.2) then rhoW = 0.d0 do k = 1, Nsp rhoW = rhoW + q(k,i)/Wk(k) enddo T = q(Nsp+3,i) qf(i) = rhoW*RU*T c # Convert dynes/cm**2 into Pa = J/m**2 if (ckunits) qf(i) = qf(i)*1.d-1 c Compute gamma else if (mcomp.eq.3) then rhoW = 0.d0 do k = 1, Nsp rhoW = rhoW + q(k,i)/Wk(k) enddo T = q(Nsp+3,i) rhoCp = avgtabip(T,q(1,i),cpk,Nsp) qf(i) = RU/(rhoCp/rhoW-RU)+1.d0 c Compute Yk else if (mcomp.ge.4.and.mcomp.lt.4+Nsp) then rho = 0.d0 do k = 1, Nsp rho = rho + q(k,i) enddo qf(i) = q(mcomp-3,i)/rho else qf(i) = 0.d0 endif 10 continue return end