subroutine remap( nleft, nright, ngeom ) C C Remap mass, momentum, and energy from the updated lagrangian grid C to the fixed Eulerian grid, using piecewise parabolic functions. C No flattening is used on remap (pass dummy array: dum to parabola.f). C include 'global.h' include 'sweep.h' real du(maxsweep), ul(maxsweep), u6(maxsweep) & ,dv(maxsweep), vl(maxsweep), v6(maxsweep) & ,dw(maxsweep), wl(maxsweep), w6(maxsweep) & ,de(maxsweep), el(maxsweep), e6(maxsweep) & ,dp(maxsweep), pl(maxsweep), p6(maxsweep) & ,dr(maxsweep), rl(maxsweep), r6(maxsweep) & ,dmu(maxsweep), dmu0(maxsweep) & ,paracoff( 9, maxsweep ) & ,delta(maxsweep), dum(maxsweep) real fluxr(maxsweep), fluxu(maxsweep), & fluxe(maxsweep), fluxv(maxsweep), fluxw(maxsweep) real rr, ur, vr, wr, pr, er, twothd, third, fractn, fractn2 integer n, ngeom C rndoff = 1.e-09 third = 1. / 3. twothd = 2. * third fourthd= 4. * third C C--------------------------------------------------------------------------- C C First calculate the volume of the overlapping subshells (delta) C and the total mass (dmu) depending on geometry C if(ngeom.eq.0) then do 1 n = nmin, nmax+1 delta(n) = xa(n) - xa0(n) 1 continue else if(ngeom.eq.1) then do 2 n = nmin, nmax+1 delta(n) = xa(n) - xa0(n) delta(n) = delta(n)*(xa0(n)+.5*delta(n)) 2 continue else if(ngeom.eq.2) then do 3 n = nmin, nmax+1 delta(n) = xa(n) - xa0(n) delta(n) = delta(n)*(xa0(n)*(xa0(n) & + delta(n)) + delta(n)**2*third) 3 continue else if(ngeom.eq.3) then do 4 n = nmin, nmax+1 delta(n) = (xa(n) - xa0(n)) * radius 4 continue else if(ngeom.eq.4) then do 5 n = nmin, nmax+1 delta(n) = (cos(xa0(n)) - cos(xa(n))) * radius 5 continue else if(ngeom.eq.5) then do 6 n = nmin, nmax+1 delta(n) = (xa(n) - xa0(n)) * radius * stheta 6 continue endif C C Generate interpolation functions, saving da, al for C constructing left and right total energy states. C call paraset( paracoff, dx, xa, nmin, nmax, ngeom ) call parabola(nmin,nmax,paracoff,rho,dr,r6,rl,dum,0,0,ngeom) call parabola(nmin,nmax,paracoff,u ,du,u6,ul,dum,0,0,ngeom) call parabola(nmin,nmax,paracoff,v ,dv,v6,vl,dum,0,0,ngeom) call parabola(nmin,nmax,paracoff,w ,dw,w6,wl,dum,0,0,ngeom) call parabola(nmin,nmax,paracoff,p ,dp,p6,pl,dum,0,0,ngeom) C C Use the profiles for density, pressure, and velocities to C calculate consistent values of the left and right values of total energy C do 12 n = nmin, nmax rr = rl(n) + dr(n) ur = ul(n) + du(n) pr = pl(n) + dp(n) vr = vl(n) + dv(n) wr = wl(n) + dw(n) el(n) = pl(n)/(rl(n)*gamm)+.5*(ul(n)**2+vl(n)**2+wl(n)**2) er = pr /(rr *gamm)+.5*(ur**2 +vr**2 +wr**2) de(n) = er-el(n) 12 continue C C Use the constructed left and right states to get parabolas for total energy C call parabola(nmin,nmax,paracoff,enrg,de,e6,el,dum,0,1,ngeom) C C Apply boundary conditions for flow onto the grid C if(nleft.eq.0 .or. nleft.eq.1) then rl(nmin-1) = rl(nmin) dr(nmin-1) = 0. r6(nmin-1) = 0. vl(nmin-1) = vl(nmin) dv(nmin-1) = 0. v6(nmin-1) = 0. wl(nmin-1) = wl(nmin) dw(nmin-1) = 0. w6(nmin-1) = 0. ul(nmin-1) = ul(nmin) du(nmin-1) = 0. u6(nmin-1) = 0. el(nmin-1) = el(nmin) de(nmin-1) = 0. e6(nmin-1) = 0. else if(nleft .eq. 2) then rl(nmin-1) = dinflo dr(nmin-1) = 0.0 r6(nmin-1) = 0.0 vl(nmin-1) = vinflo dv(nmin-1) = 0.0 v6(nmin-1) = 0.0 wl(nmin-1) = winflo dw(nmin-1) = 0.0 w6(nmin-1) = 0.0 ul(nmin-1) = uinflo du(nmin-1) = 0.0 u6(nmin-1) = 0.0 el(nmin-1) = einflo de(nmin-1) = 0.0 e6(nmin-1) = 0.0 else if(nleft .eq. 3) then rl(nmin-1) = rl(nmax) dr(nmin-1) = dr(nmax) r6(nmin-1) = r6(nmax) vl(nmin-1) = vl(nmax) dv(nmin-1) = dv(nmax) v6(nmin-1) = v6(nmax) wl(nmin-1) = wl(nmax) dw(nmin-1) = dw(nmax) w6(nmin-1) = w6(nmax) ul(nmin-1) = ul(nmax) du(nmin-1) = du(nmax) u6(nmin-1) = u6(nmax) el(nmin-1) = el(nmax) de(nmin-1) = de(nmax) e6(nmin-1) = e6(nmax) else if(nleft .eq. 5) then rl(nmin-1) = dinflo dr(nmin-1) = 0.0 r6(nmin-1) = 0.0 vl(nmin-1) = winflo dv(nmin-1) = 0.0 v6(nmin-1) = 0.0 wl(nmin-1) = uinflo dw(nmin-1) = 0.0 w6(nmin-1) = 0.0 ul(nmin-1) = vinflo du(nmin-1) = 0.0 u6(nmin-1) = 0.0 el(nmin-1) = einflo de(nmin-1) = 0.0 e6(nmin-1) = 0.0 endif if(nright.eq.0 .or. nright.eq.1) then rl(nmax+1) = rl(nmax) dr(nmax+1) = 0. r6(nmax+1) = 0. vl(nmax+1) = vl(nmax) dv(nmax+1) = 0. v6(nmax+1) = 0. wl(nmax+1) = wl(nmax) dw(nmax+1) = 0. w6(nmax+1) = 0. ul(nmax+1) = ul(nmax) du(nmax+1) = 0. u6(nmax+1) = 0. el(nmax+1) = el(nmax) de(nmax+1) = 0. e6(nmax+1) = 0. else if(nright .eq. 2) then rl(nmax+1) = dinflo dr(nmax+1) = 0.0 r6(nmax+1) = 0.0 vl(nmax+1) = vinflo dv(nmax+1) = 0.0 v6(nmax+1) = 0.0 wl(nmax+1) = winflo dw(nmax+1) = 0.0 w6(nmax+1) = 0.0 ul(nmax+1) = uinflo du(nmax+1) = 0.0 u6(nmax+1) = 0.0 el(nmax+1) = einflo de(nmax+1) = 0.0 e6(nmax+1) = 0.0 else if(nright .eq. 3) then rl(nmax+1) = rl(nmin) dr(nmax+1) = dr(nmin) r6(nmax+1) = r6(nmin) vl(nmax+1) = vl(nmin) dv(nmax+1) = dv(nmin) v6(nmax+1) = v6(nmin) wl(nmax+1) = wl(nmin) dw(nmax+1) = dw(nmin) w6(nmax+1) = w6(nmin) ul(nmax+1) = ul(nmin) du(nmax+1) = du(nmin) u6(nmax+1) = u6(nmin) el(nmax+1) = el(nmin) de(nmax+1) = de(nmin) e6(nmax+1) = e6(nmin) else if(nright .eq. 5) then rl(nmax+1) = dinflo dr(nmax+1) = 0. r6(nmax+1) = 0. vl(nmax+1) = winflo dv(nmax+1) = 0. v6(nmax+1) = 0. wl(nmax+1) = uinflo dw(nmax+1) = 0. w6(nmax+1) = 0. ul(nmax+1) = vinflo du(nmax+1) = 0. u6(nmax+1) = 0. el(nmax+1) = einflo de(nmax+1) = 0. e6(nmax+1) = 0. endif C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C C Calculate the total mass (fluxr), momentum (fluxu), and energy (fluxe) C in the subshell created by the overlap of the Lagrangian and Eulerican grids. C If the zone face has moved to the left (deltx > 0), use the integral from the C left side of zone n (fluxrr). If the zone face has moved to the right C (deltx < 0), use the integral from the right side of zone k=n-1 (fluxrl). C do 100 n = nmin, nmax + 1 deltx = xa(n) - xa0(n) if(deltx .ge. 0.0) then k = n - 1 fractn = 0.5*deltx/dx(k) fractn2 = 1. - fourthd*fractn fluxr(n) = delta(n)*( & rl(k)+dr(k)-fractn*(dr(k)-fractn2*r6(k))) fluxu(n) = ul(k)+du(k)-fractn*(du(k)-fractn2*u6(k)) fluxv(n) = vl(k)+dv(k)-fractn*(dv(k)-fractn2*v6(k)) fluxw(n) = wl(k)+dw(k)-fractn*(dw(k)-fractn2*w6(k)) fluxe(n) = el(k)+de(k)-fractn*(de(k)-fractn2*e6(k)) else fractn = 0.5*deltx/dx(n) fractn2 = 1. + fourthd*fractn fluxr(n) = delta(n)*( & rl(n) - fractn*(dr(n)+fractn2*r6(n))) fluxu(n) = ul(n) - fractn*(du(n)+fractn2*u6(n)) fluxv(n) = vl(n) - fractn*(dv(n)+fractn2*v6(n)) fluxw(n) = wl(n) - fractn*(dw(n)+fractn2*w6(n)) fluxe(n) = el(n) - fractn*(de(n)+fractn2*e6(n)) endif 100 continue C C Advect mass, momentum, and energy by moving the subshell quantities C into the appropriate Eulerian zone. C do 120 n = nmin, nmax dmu (n) = rho(n) * dvol(n) dmu0(n) = (dmu(n) + fluxr(n) - fluxr(n+1)) rho (n) = dmu0(n)/dvol0(n) rho (n) = max(smallr,rho(n)) dmu0(n) = 1./dmu0(n) c dlogdm = abs(dmu(n)-1./dmu0(n))*dmu0(n) c if(dlogdm.lt.rndoff) dmu0(n) = dmu(n) u (n) = (u (n)*dmu(n) + (fluxr(n)*fluxu(n) & - fluxr(n+1)*fluxu(n+1)) )*dmu0(n) v (n) = (v (n)*dmu(n) + (fluxr(n)*fluxv(n) & - fluxr(n+1)*fluxv(n+1)) )*dmu0(n) w (n) = (w (n)*dmu(n) + (fluxr(n)*fluxw(n) & - fluxr(n+1)*fluxw(n+1)) )*dmu0(n) enrg(n) = (enrg(n)*dmu(n) & + fluxr(n)*fluxe(n) - fluxr(n+1)*fluxe(n+1) & )*dmu0(n) c pold = p(n) p(n) = gamm*rho(n)*(enrg(n)-0.5*(u(n)**2+v(n)**2+w(n)**2)) c dlogp = abs(p(n)-pold)/pold c if(dlogp.lt.rndoff) p(n) = pold(n) p(n) = max(smallp,p(n)) 120 continue C return end