c c ===================================================== subroutine setaux(maxmx,maxmy,maxmz,meqn,mbc,ibx,iby,ibz, & mx,my,mz,q,aux,maux,cornx,corny,cornz,dx,dy,dz,t,dt) c ===================================================== c c # alter total energy and density to ensure admissible c # vectors of state c c Copyright (C) 2008 California Institute of Technology c Ralf Deiterding, ralf@cacr.caltech.edu c implicit double precision (a-h, o-z) dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc, & 1-iby*mbc:maxmy+iby*mbc, 1-ibz*mbc:maxmz+ibz*mbc) dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc, & 1-iby*mbc:maxmy+iby*mbc, 1-ibz*mbc:maxmz+ibz*mbc) c do 10 k=1-ibz*mbc,mz+ibz*mbc do 10 j=1-iby*mbc,my+iby*mbc do 10 i=1-ibx*mbc,mx+ibx*mbc if (q(1,i,j,k).lt.1.d-8) q(1,i,j,k)=1.d-8 rho = q(1,i,j,k) u = q(2,i,j,k)/q(1,i,j,k) v = q(3,i,j,k)/q(1,i,j,k) w = q(4,i,j,k)/q(1,i,j,k) p = (q(5,i,j,k) - 0.5d0*rho*(u**2+v**2+w**2) - & q(7,i,j,k))/q(6,i,j,k) gamma1 = 1.d0/q(6,i,j,k) gamma = gamma1 + 1.d0 pinfty = q(7,i,j,k)*gamma1/gamma plimit = -pinfty + pinfty*1.d-8 if (p.lt.plimit) q(5,i,j,k) = plimit*q(6,i,j,k)+ & 0.5d0*rho*(u**2+v**2+w**2)+q(7,i,j,k) 10 continue c return end