c c c ================================================================== subroutine rpt3meu(ixyz,icoor,maxm,meqn,mwaves,mbc,mx, & ql,qr,maux,aux1,aux2,aux3,ilr,asdq, & bmasdq,bpasdq) c ================================================================== implicit double precision(a-h,o-z) c c # Riemann solver in the transverse direction for the c # Euler equations. c # c # c # On input, c c # ql,qr is the data along some one-dimensional slice, as in rpn3 c # This slice is c # in the x-direction if ixyz=1, c # in the y-direction if ixyz=2, or c # in the z-direction if ixyz=3. c # asdq is an array of flux differences (A^* \Delta q). c # asdq(i,:) is the flux difference propagating away from c # the interface between cells i-1 and i. c # imp = 1 if asdq = A^- \Delta q, the left-going flux difference c # 2 if asdq = A^+ \Delta q, the right-going flux difference c c # aux2 is the auxiliary array (if method(7)=maux>0) along c # the plane where this slice lies, say at j=J if ixyz=1. c # aux2(:,:,1) contains data along j=J, k=k-1 c # aux2(:,:,2) contains data along j=J, k=k c # aux2(:,:,3) contains data along j=J, k=k+1 c # aux1 is the auxiliary array along the plane with j=J-1 c # aux3 is the auxiliary array along the plane with j=J+1 c c # if ixyz=2 then aux2 is in some plane k=K, and c # aux2(:,:,1) contains data along i=I-1, k=K, etc. c c # if ixyz=3 then aux2 is in some plane i=I, and c # aux2(:,:,1) contains data along j=j-1, i=I, etc. c c # On output, c # If data is in x-direction (ixyz=1) then this routine does the c # splitting of asdq (= A^* \Delta q, where * = + or -) ... c c # into down-going flux difference bmasdq (= B^- A^* \Delta q) c # and up-going flux difference bpasdq (= B^+ A^* \Delta q) c # when icoor = 2, c c # or c c # into down-going flux difference bmasdq (= C^- A^* \Delta q) c # and up-going flux difference bpasdq (= C^+ A^* \Delta q) c # when icoor = 3. c # c c # More generally, ixyz specifies what direction the slice of data is c # in, and icoor tells which transverse direction to do the splitting in: c c # If ixyz = 1, data is in x-direction and then c # icoor = 2 => split in the y-direction c # icoor = 3 => split in the z-direction c c # If ixyz = 2, data is in y-direction and then c # icoor = 2 => split in the z-direction c # icoor = 3 => split in the x-direction c c # If ixyz = 3, data is in z-direction and then c # icoor = 2 => split in the x-direction c # icoor = 3 => split in the y-direction c c # c # Uses Roe averages and other quantities which were c # computed in rpn3eu and stored in the common block comroe. c c # Copyright (C) 2002 Ralf Deiterding c # Brandenburgische Universitaet Cottbus c dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension asdq(1-mbc:maxm+mbc, meqn) dimension bmasdq(1-mbc:maxm+mbc, meqn) dimension bpasdq(1-mbc:maxm+mbc, meqn) dimension aux1(1-mbc:maxm+mbc, maux, 3) dimension aux2(1-mbc:maxm+mbc, maux, 3) dimension aux3(1-mbc:maxm+mbc, maux, 3) c dimension waveb(7,3),sb(3) parameter (maxmrp = 1005) !# assumes atmost max(mx,my,mz) = 1000 with mbc=5 parameter (minmrp = -4) !# assumes at most mbc=5 common /comroe/ u2v2w2(minmrp:maxmrp), & u(minmrp:maxmrp),v(minmrp:maxmrp),w(minmrp:maxmrp), & enth(minmrp:maxmrp),a(minmrp:maxmrp),g1a2(minmrp:maxmrp), & euv(minmrp:maxmrp),p(minmrp:maxmrp) c if (minmrp.gt.1-mbc .or. maxmrp .lt. maxm+mbc) then write(6,*) 'need to increase maxmrp in rp3t' stop endif c if(ixyz .eq. 1)then mu = 2 mv = 3 mw = 4 else if(ixyz .eq. 2)then mu = 3 mv = 4 mw = 2 else mu = 4 mv = 2 mw = 3 endif c c # Solve Riemann problem in the second coordinate direction c if( icoor .eq. 2 )then do 20 i = 2-mbc, mx+mbc a4 = g1a2(i) * (euv(i)*asdq(i,1) & + u(i)*asdq(i,mu) + v(i)*asdq(i,mv) & + w(i)*asdq(i,mw) - asdq(i,5) & + p(i)*asdq(i,6) + asdq(i,7) ) a2 = asdq(i,mu) - u(i)*asdq(i,1) a3 = asdq(i,mw) - w(i)*asdq(i,1) a5 = (asdq(i,mv) + (a(i)-v(i))*asdq(i,1) - a(i)*a4) & / (2.d0*a(i)) a1 = asdq(i,1) - a4 - a5 a6 = asdq(i,6) a7 = asdq(i,7) c waveb(1,1) = a1 waveb(mu,1) = a1*u(i) waveb(mv,1) = a1*(v(i)-a(i)) waveb(mw,1) = a1*w(i) waveb(5,1) = a1*(enth(i) - v(i)*a(i)) waveb(6,1) = 0.d0 waveb(7,1) = 0.d0 sb(1) = v(i) - a(i) c waveb(1,2) = a4 waveb(mu,2) = a2 + u(i)*a4 waveb(mv,2) = v(i)*a4 waveb(mw,2) = a3 + w(i)*a4 waveb(5,2) = a4*0.5d0*u2v2w2(i)+a2*u(i)+a3*w(i)+a6*p(i)+a7 waveb(6,2) = a6 waveb(7,2) = a7 sb(2) = v(i) c waveb(1,3) = a5 waveb(mu,3) = a5*u(i) waveb(mv,3) = a5*(v(i)+a(i)) waveb(mw,3) = a5*w(i) waveb(5,3) = a5*(enth(i)+v(i)*a(i)) waveb(6,3) = 0.d0 waveb(7,3) = 0.d0 sb(3) = v(i) + a(i) c do 25 m=1,meqn bmasdq(i,m) = 0.d0 bpasdq(i,m) = 0.d0 do 25 mws=1,mwaves bmasdq(i,m) = bmasdq(i,m) & + dmin1(sb(mws), 0.d0) * waveb(m,mws) bpasdq(i,m) = bpasdq(i,m) & + dmax1(sb(mws), 0.d0) * waveb(m,mws) 25 continue c 20 continue c else c c # Solve Riemann problem in the third coordinate direction c do 30 i = 2-mbc, mx+mbc a4 = g1a2(i) * (euv(i)*asdq(i,1) & + u(i)*asdq(i,mu) + v(i)*asdq(i,mv) & + w(i)*asdq(i,mw) - asdq(i,5) & + p(i)*asdq(i,6) + asdq(i,7) ) a2 = asdq(i,mu) - u(i)*asdq(i,1) a3 = asdq(i,mv) - v(i)*asdq(i,1) a5 = (asdq(i,mw) + (a(i)-w(i))*asdq(i,1) - a(i)*a4) & / (2.d0*a(i)) a1 = asdq(i,1) - a4 - a5 a6 = asdq(i,6) a7 = asdq(i,7) c waveb(1,1) = a1 waveb(mu,1) = a1*u(i) waveb(mv,1) = a1*v(i) waveb(mw,1) = a1*(w(i) - a(i)) waveb(5,1) = a1*(enth(i) - w(i)*a(i)) waveb(6,1) = 0.d0 waveb(7,1) = 0.d0 sb(1) = w(i) - a(i) c waveb(1,2) = a4 waveb(mu,2) = a2 + u(i)*a4 waveb(mv,2) = a3 + v(i)*a4 waveb(mw,2) = w(i)*a4 waveb(5,2) = a4*0.5d0*u2v2w2(i)+a2*u(i)+a3*v(i)+a6*p(i)+a7 waveb(6,2) = a6 waveb(7,2) = a7 sb(2) = w(i) c waveb(1,3) = a5 waveb(mu,3) = a5*u(i) waveb(mv,3) = a5*v(i) waveb(mw,3) = a5*(w(i)+a(i)) waveb(5,3) = a5*(enth(i)+w(i)*a(i)) waveb(6,3) = 0.d0 waveb(7,3) = 0.d0 sb(3) = w(i) + a(i) c do 35 m=1,meqn bmasdq(i,m) = 0.d0 bpasdq(i,m) = 0.d0 do 35 mws=1,mwaves bmasdq(i,m) = bmasdq(i,m) & + dmin1(sb(mws), 0.d0) * waveb(m,mws) bpasdq(i,m) = bpasdq(i,m) & + dmax1(sb(mws), 0.d0) * waveb(m,mws) 35 continue c 30 continue c endif c return end