Describe the proposal
Add the face-centered velocities returned by the Riemann solver to the global state MultiFabs state_fc_new_ and state_fc_old_. (Note that these velocities do not have the same time centering as other variables. They are always at the midpoint of the timestep, since they reflect the average face-centered velocity over the timestep.)

This is necessary to implement tracer particles that are advected with the same velocities as the fluid. This requires using the face-centered velocities from the Riemann solver. Ghost cells (ghost faces?) that must be filled with AMR interpolation are required because particles may cross the coarse-fine boundary during a timestep.

Caveat: This causes a slight complication for users, because simulations must now initialize the face-centered velocity field as part of the initial conditions or the code will complain (or possibly crash). If tracer particles are used, boundary conditions for the face-centered velocity must be provided (or the code will crash).

Describe alternatives you've considered
The alternative is to use cell-centered velocities for tracer particle advection. However, this is not exactly consistent with the mass flux of the fluid. While this difference may seem small, it can lead to order-of-magnitude errors in tracer particle advection in supersonic turbulence (see https://ui.adsabs.harvard.edu/abs/2013MNRAS.435.1426G/abstract and http://girichidis.de/ringberg2016/Genel.pdf).

Additional context
This is implemented in #490.