                     dotr += wave[mq]*waves[I_waves+1];
                 }
                 I_speeds = I + mw*zs;
-                r = (speeds[I_speeds] > 0) ? dotl/wnorm2 : dotr/wnorm2;
-
-                wlimitr = limiter(r);  /* allow for selection */
-
+                s[mw] = speeds[I_speeds];
+                wlimitr = 1;
+                if (wnorm2 != 0)
+                {
+                    r = (s[mw] > 0) ? dotl/wnorm2 : dotr/wnorm2;
+                    wlimitr = limiter(r);  /* allow for selection */
+                }
+
+#if 0
                 for(mq = 0; mq < meqn; mq++)
                 {
                     I_waves = I + (mw*meqn + mq)*zs;
                     waves[I_waves] = wlimitr*wave[mq];                    
                 }
+#endif                
+                for(mq = 0; mq < meqn; mq++)
+                {
+                    I_q = I + mq;
+                    cqxx += fabs(s[mw])*(1.0 - fabs(s[mw])*dtdx)*wlimitr*wave[mq];
+                    fm[I_q] += 0.5*cqxx;   /* amdq + cqxx */
+                    fp[I_q] += 0.5*cqxx;   /* apdq - cqxx */
+                }
+
 
                 /* y-faces */
                 wnorm2 = 0;

        I = (ix + mbc)*xs + (iy + mbc)*ys;

        double wnorm2,dotr,dotl, wlimitr,r;
        double wave[20];
        if (ix < mx + 1 && iy < my + 1)
        {
            for(mw = 0; mw < mwaves; mw++)
            {
                /* x-faces */

                    dotr += wave[mq]*waves[I_waves+ys];
                }
                I_speeds = I + (mwaves + mw)*zs;
                s[mw] = speeds[I_speeds];
                wlimitr = 1;
                if (wnorm2 != 0)
                {
                    r = (s[mw] > 0) ? dotl/wnorm2 : dotr/wnorm2;
                    wlimitr = limiter(r);  /* allow for selection */
                }

#if 0
                for(mq = 0; mq < meqn; mq++)
                {
                    I_waves = I + ((mwaves+mw)*meqn + mq)*zs;
                    waves[I_waves] = wlimitr*wave[mq];                    
                }    
#endif
                for(mq = 0; mq < meqn; mq++)
                {
                    I_q = I + mq;
                    cqxx += fabs(s[mw])*(1.0 - fabs(s[mw])*dtdx)*wlimitr*wave[mq];
                    gm[I_q] += 0.5*cqxx;   /* amdq + cqxx */
                    gp[I_q] += 0.5*cqxx;   /* apdq - cqxx */
                }


            }
        }
    }
