c c Copyright (C) 2002 Ralf Deiterding c Brandenburgische Universitaet Cottbus c c ========================================================= subroutine src(maxmx,meqn,mbc,ibx,mx,q, & aux,maux,t,dt,ibnd) c ========================================================= implicit double precision(a-h,o-z) include "cuser.i" c dimension q(meqn, 1-ibx*mbc:maxmx+ibx*mbc) dimension aux(maux, 1-ibx*mbc:maxmx+ibx*mbc) c dimension y(1),c(24),w(1,9) common /euparam/ e external fcn3 c n=1 tol=1.d-12 nw=1 pmax = 0.d0 do 10 i=1-ibnd*ibx*mbc,mx+ibnd*ibx*mbc rho=q(1,i)+q(2,i) u=q(3,i)/rho e=q(4,i)/rho-0.5d0*u**2 temp=gamma1*(e-q(2,i)/rho*q0) if (temp.gt.1.d0) then ind=1 z=t zend=t+dt y(1)=q(2,i)/rho call dverk(n,fcn3,z,y,zend,tol,ind,c,nw,w) c do n=1,100 c call fcn3(1,x,y,yprime) c y(1) = y(1)+dt*yprime/100.d0 c enddo q(2,i)=rho*y(1) q(1,i)=rho-q(2,i) endif 10 continue c return end c c ========================================================= subroutine fcn3(n,x,y,yprime) c ========================================================= implicit double precision (a-h,o-z) include "cuser.i" c dimension y(n),yprime(n) common /euparam/ e c temp=gamma1*(e-y(1)*q0) yprime(1)=-dhalbe*y(1)*exp(-E0/temp) return end