diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py
index 073d86624..d81501263 100644
--- a/pyro/compressible_fv4/fluxes.py
+++ b/pyro/compressible_fv4/fluxes.py
@@ -50,7 +50,7 @@ def fluxes(myd, rp, ivars):
     # get the cell-average data
     U_avg = myd.data
 
-    # convert U from cell-centers to cell averages
+    # convert U from cell-averages to cell-centers
     U_cc = np.zeros_like(U_avg)
 
     U_cc[:, :, ivars.idens] = myd.to_centers("density")
@@ -81,6 +81,8 @@ def fluxes(myd, rp, ivars):
     # for debugging
     nolimit = 0
 
+    # do reconstruction on the cell-average primitive variable state, q_avg
+
     for idir in [1, 2]:
 
         # interpolate <W> to faces (with limiting)
@@ -117,6 +119,7 @@ def fluxes(myd, rp, ivars):
                     q_r.v(n=n, buf=2)[:, :] = xi.v(buf=2)*q_r.v(n=n, buf=2) + \
                         (1.0 - xi.v(buf=2))*q_avg.v(n=n, buf=2)
 
+        # solve the Riemann problem to find the face-average q
         _q = cf.riemann_prim(idir, myg.ng,
                              ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux,
                              0, 0,
@@ -124,7 +127,7 @@ def fluxes(myd, rp, ivars):
 
         q_int_avg = ai.ArrayIndexer(_q, grid=myg)
 
-        # calculate the face-centered W using the transverse Laplacian
+        # calculate the face-centered q using the transverse Laplacian
         q_int_fc = myg.scratch_array(nvar=ivars.nq)
 
         if idir == 1:
@@ -140,7 +143,8 @@ def fluxes(myd, rp, ivars):
                                 2*q_int_avg.v(n=n, buf=myg.ng-1) +
                                 q_int_avg.ip(-1, n=n, buf=myg.ng-1))
 
-        # compute the final fluxes
+        # compute the final fluxes using both the face-average state, q_int_avg,
+        # and face-centered q, q_int_fc
         F_fc = flux_cons(ivars, idir, gamma, q_int_fc)
         F_avg = flux_cons(ivars, idir, gamma, q_int_avg)