                         wave[mq] = waves[I_waves];
                     }                        
 
-                    //if (mthlim[mw] > 0)
+                    if (mthlim[mw] > 0)
                     {
                         wnorm2 = dotl = dotr = 0;
                         for(mq = 0; mq < meqn; mq++)
                         {
                             I_waves = I + ((mwaves+mw)*meqn + mq)*zs;
                             wnorm2 += pow(wave[mq],2);
-                            dotl += wave[mq]*waves[I_waves-ys];
-                            dotr += wave[mq]*waves[I_waves+ys];
-                        }
+                        }  
                         if (wnorm2 != 0)
                         {
+                            for(mq = 0; mq < meqn; mq++)
+                            {
+                                I_waves = I + ((mwaves+mw)*meqn + mq)*zs;
+                                dotl += wave[mq]*waves[I_waves-ys];
+                                dotr += wave[mq]*waves[I_waves+ys];
+                            }  
                             r = (s[mw] > 0) ? dotl/wnorm2 : dotr/wnorm2;
+
                             wlimitr = cudaclaw_limiter(mthlim[mw],r);  
-                        }
-                        for (mq = 0; mq < meqn; mq++)
-                        {
-                            wave[mq] *= wlimitr;
+
+                            for (mq = 0; mq < meqn; mq++)
+                            {
+                                wave[mq] *= wlimitr;
+                            }
                         }
                     }
 

                    {
                        I_q = I + mq*zs;
                        cqyy = fabs(s[mw])*(1.0 - fabs(s[mw])*dtdy)*wave[mq];
                        gm[I_q] += 0.5*cqyy;   
                        gp[I_q] += 0.5*cqyy;  
                        //if (order[1] > 0)
                        {                            
                            bmdq_trans[I_q] += cqyy;     
                            bpdq_trans[I_q] -= cqyy;      

        ix = thread_index % ifaces_x;
        iy = thread_index/ifaces_y;

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

        //if (ix < mx && 0 < iy && iy < my+1)   /* Is this needed? */
        {
            for(mq = 0; mq < meqn; mq++)
            {
