c c c ===================================================== subroutine rpn3euznd(ixyz,maxm,meqn,mwaves,mbc,mx,ql,qr, & maux,auxl,auxr,wave,s,fl,fr) c ===================================================== c c # solve Riemann problems for the 3D ZND-Euler equations using c # van Leer's Flux Vector Splitting following Shuen's approach for c # multicomponent gas mixtures c c # On input, ql contains the state vector at the left edge of each cell c # qr contains the state vector at the right edge of each cell c # This data is along a slice in the x-direction if ixyz=1 c # the y-direction if ixyz=2. c # the z-direction if ixyz=3. c c # On output, wave contains the waves, s the speeds, c # fl and fr the positive and negative flux. c c # Note that the i'th Riemann problem has left state qr(i-1,:) c # and right state ql(i,:) c # From the basic clawpack routines, this routine is called with ql = qr c c # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c implicit double precision (a-h,o-z) c dimension wave(1-mbc:maxm+mbc, meqn, mwaves) dimension s(1-mbc:maxm+mbc, mwaves) dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension fl(1-mbc:maxm+mbc, meqn) dimension fr(1-mbc:maxm+mbc, meqn) dimension auxl(1-mbc:maxm+mbc, maux, 3) dimension auxr(1-mbc:maxm+mbc, maux, 3) double precision Ml, Mr dimension fvl(6), fvr(6), sl(3), sr(3) common /param/ gamma,gamma1,q0 c c # Method returns fluxes c ------------ common /rpnflx/ mrpnflx mrpnflx = 1 c c # set mu to point to the component of the system that corresponds c # to momentum in the direction of this slice, mv and mw to the c # orthogonal momentums: c if(ixyz .eq. 1)then mu = 3 mv = 4 mw = 5 else if(ixyz .eq. 2)then mu = 4 mv = 5 mw = 3 else mu = 5 mv = 3 mw = 4 endif c c # Van Leer's Flux Vector Splitting c gamma2 = gamma**2-1 do 10 i=2-mbc,mx+mbc rhol = qr(i-1,1)+qr(i-1,2) rhor = ql(i ,1)+ql(i ,2) Y1l = qr(i-1,1)/rhol Y2l = qr(i-1,2)/rhol Y1r = ql(i ,1)/rhor Y2r = ql(i ,2)/rhor ul = qr(i-1,mu)/rhol ur = ql(i ,mu)/rhor vl = qr(i-1,mv)/rhol vr = ql(i ,mv)/rhor wl = qr(i-1,mw)/rhol wr = ql(i ,mw)/rhor El = qr(i-1,6)/rhol Er = ql(i ,6)/rhor pl = gamma1*(qr(i-1,6) - qr(i-1,2)*q0 - & 0.5d0*(qr(i-1,mu)**2+qr(i-1,mv)**2+qr(i-1,mw)**2)/rhol) pr = gamma1*(ql(i ,6) - ql(i ,2)*q0 - & 0.5d0*(ql(i ,mu)**2+ql(i ,mv)**2+ql(i ,mw)**2)/rhor) Hl = El+pl/rhol Hr = Er+pr/rhor c al2 = gamma*pl/rhol al = dsqrt(al2) ar2 = gamma*pr/rhor ar = dsqrt(ar2) c Ml = ul/al Mr = ur/ar c sl(1) = ul-al sl(2) = ul sl(3) = ul+al sr(1) = ur-ar sr(2) = ur sr(3) = ur+ar c if (Ml.ge.1.d0) then facl = rhol*ul fvl(1) = facl*Y1l fvl(2) = facl*Y2l fvl(mu) = facl*ul+pl fvl(mv) = facl*vl fvl(mw) = facl*wl fvl(6) = facl*El+ul*pl else if (Ml.le.-1.d0) then do m = 1,meqn fvl(m) = 0.d0 enddo else fhl = gamma*(El - Y2l*q0 - 0.5d0*(ul**2+vl**2+wl**2))/al2 xl = fhl/(1.d0+2.d0*fhl) facl = 0.25d0*rhol*al*(Ml+1.d0)**2 fvl(1) = facl*Y1l fvl(2) = facl*Y2l fvl(mu) = facl*2.d0*al/gamma*(0.5d0*gamma1*Ml+1.d0) fvl(mv) = facl*vl fvl(mw) = facl*wl fvl(6) = facl*(Hl-xl*(ul-al)**2) endif c if (Mr.le.-1.d0) then facr = rhor*ur fvr(1) = facr*Y1r fvr(2) = facr*Y2r fvr(mu) = facr*ur+pr fvr(mv) = facr*vr fvr(mw) = facr*wr fvr(6) = facr*Er+ur*pr else if (Mr.ge.1.d0) then do m = 1,meqn fvr(m) = 0.d0 enddo else fhr = gamma*(Er - Y2r*q0 - 0.5d0*(ur**2+vr**2+wr**2))/ar2 xr = fhr/(1.d0+2.d0*fhr) facr = -0.25d0*rhor*ar*(Mr-1.d0)**2 fvr(1) = facr*Y1r fvr(2) = facr*Y2r fvr(mu) = facr*2.d0*ar/gamma*(0.5d0*gamma1*Mr-1.d0) fvr(mv) = facr*vr fvr(mw) = facr*wr fvr(6) = facr*(Hr-xr*(ur+ar)**2) endif c do 20 m = 1,meqn fl(i,m) = fvl(m) + fvr(m) fr(i,m) = -fl(i,m) 20 continue c if (dabs(Ml).lt.1.d0) then facl = (gamma+3.d0)/(2.d0*gamma+dabs(Ml)*(3.d0-gamma)) else facl = 1.d0 endif if (dabs(Mr).lt.1.d0) then facr = (gamma+3.d0)/(2.d0*gamma+dabs(Mr)*(3.d0-gamma)) else facr = 1.d0 endif c do 10 mws=1,mwaves s(i,mws) = dmax1(dabs(facl*sl(mws)),dabs(facr*sr(mws))) do 10 m=1,meqn wave(i,m,mws) = 0.d0 10 continue c return end c