         ifaces_y = my + 1;
         num_ifaces = ifaces_x*ifaces_y;
 
-        for(thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
+        for(int thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
         { 
-            ix = thread_index % ifaces_x;
-            iy = thread_index/ifaces_x;
+            int ix = thread_index % ifaces_x;
+            int iy = thread_index/ifaces_x;
 
             /* Start at first non-ghost interior cell */
-            I = (iy + mbc)*ys + ix + mbc - 1;
-
-            for(mw = 0; mw < mwaves; mw++)
+            int I = (iy + mbc)*ys + ix + mbc - 1;
+
+            for(int mw = 0; mw < mwaves; mw++)
             {
                 /* ------------------------------- Y-directions --------------------------- */
-                I_speeds = I + (mwaves + mw)*zs;
+                int I_speeds = I + (mwaves + mw)*zs;
                 s[mw] = speeds[I_speeds];
 
-                for(mq = 0; mq < meqn; mq++)
+                for(int mq = 0; mq < meqn; mq++)
                 {
-                    I_waves = I + ((mwaves + mw)*meqn + mq)*zs;
+                    int I_waves = I + ((mwaves + mw)*meqn + mq)*zs;
                     wave[mq] = waves[I_waves];
                 }                        
 
                 if (mthlim[mw] > 0)
                 {
-                    wnorm2 = dotl = dotr = 0;
-                    for(mq = 0; mq < meqn; mq++)
+                    double wnorm2 = 0, dotl = 0, dotr = 0;
+                    for(int mq = 0; mq < meqn; mq++)
                     {
                         wnorm2 += pow(wave[mq],2);
-                        I_waves = I + ((mwaves + mw)*meqn + mq)*zs;
+                        int I_waves = I + ((mwaves + mw)*meqn + mq)*zs;
                         dotl += wave[mq]*waves[I_waves-ys];
                         dotr += wave[mq]*waves[I_waves+ys];
                     }  
                     wnorm2 = (wnorm2 == 0) ? 1e-15 : wnorm2;  
 
-                    r = (s[mw] > 0) ? dotl/wnorm2 : dotr/wnorm2;
-
-                    wlimitr = cudaclaw_limiter(mthlim[mw],r);  
-
-                    for (mq = 0; mq < meqn; mq++)
+                    double r = (s[mw] > 0) ? dotl/wnorm2 : dotr/wnorm2;
+
+                    double wlimitr = cudaclaw_limiter(mthlim[mw],r);  
+
+                    for (int mq = 0; mq < meqn; mq++)
                     {
                         wave[mq] *= wlimitr;
                     }                    
                 }
 
-                for(mq = 0; mq < meqn; mq++)
+                for(int mq = 0; mq < meqn; mq++)
                 {
-                    I_q = I + mq*zs;
-                    if (use_fwaves != 0)
-                    {
-                        cqyy = copysign(1.,s[mw])*(1.0 - fabs(s[mw])*dtdx)*wave[mq];
-                    }
-                    else
-                    {
-                        cqyy = fabs(s[mw])*(1.0 - fabs(s[mw])*dtdx)*wave[mq];
-                    }
+                    int I_q = I + mq*zs;
+                    double cqyy = (1.0 - fabs(s[mw])*dtdx)*wave[mq];
+                    cqyy *= (use_fwaves) ? copysign(1.,s[mw]) : fabs(s[mw]);
 
                     gm[I_q] += 0.5*cqyy;   
                     gp[I_q] += 0.5*cqyy;  

            }
        }

        for(mw = 0; mw < mwaves; mw++)
        {
            maxcfl = max(maxcfl,abs(s[mw]*dtdx));

            if (order[0] == 2)
            {                    
                I_speeds = I + mw*zs;
                speeds[I_speeds] = s[mw];
                for(mq = 0; mq < meqn; mq++)
                {
                    k = mw*meqn + mq;
                    I_waves = I + k*zs;
                    waves[I_waves] = wave[k];
                }
            }


    /* ------------------------------- Final update ----------------------------------- */

    for(int thread_index = threadIdx.x; thread_index < mx*my; thread_index += blockDim.x)
    {
        int ix = thread_index % mx;
        int iy = thread_index/mx;

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

        for(int mq = 0; mq < meqn; mq++)
        {
            int I_q = I + mq*zs;
            qold[I_q] = qold[I_q] - dtdx * (fm[I_q + 1] - fp[I_q]) 
                                  - dtdy * (gm[I_q +ys] - gp[I_q]);
        }        
