     ifaces_x = mx + 1;
     ifaces_y = my + 1;
     num_ifaces = ifaces_x*ifaces_y;
- 
+
     for(int thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
     { 
         int ix = thread_index % ifaces_x;

        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+1];
                }
                I_speeds = I + 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 + (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;
