             gp[I_q - 1] -= gupdate;   
         }            
 
-        //__syncthreads();
+        __threadfence_block();
 
         for(mq = 0; mq < meqn; mq++)
         {


        I = (iy + 1)*ys + (ix + 1);  /* Start one cell from left/bottom edge */

        //if (ix < mx + 2*mbc-1 && iy < my + 2*mbc-1)
        {
            for(mq = 0; mq < meqn; mq++)
            {
                I_q = I + mq*zs;
                qr[mq] = qold[I_q];        /* Right */
                ql[mq] = qold[I_q - 1];    /* Left  */
                qd[mq] = qold[I_q - ys];   /* Down  */  
            }

            for(m = 0; m < maux; m++)
            {
                I_aux = I + m*zs;
                auxl[m] = aux[I_aux - 1];
                auxr[m] = aux[I_aux];
                auxd[m] = aux[I_aux - ys];
            }                        

            /* ------------------------ Normal solve in X direction ------------------- */
            
            rpn2(0, meqn, mwaves, maux, ql, qr, auxl, auxr, wave, s, amdq, apdq);

            for (mq = 0; mq < meqn; mq++) 
            {
                I_q = I + mq*zs;
                fm[I_q] = amdq[mq];
                fp[I_q] = -apdq[mq]; 
                if (order[1] > 0)
                {
                    amdq_trans[I_q] = amdq[mq];                                        
                    apdq_trans[I_q] = apdq[mq];  
                }
            }

            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];
                    }
                }
            }
        
            /* ------------------------ Normal solve in Y direction ------------------- */
            rpn2(1, meqn, mwaves, maux, qd, qr, auxd, auxr, wave, s, bmdq, bpdq);

            /* Set value at bottom interface of cell I */
            for (mq = 0; mq < meqn; mq++) 
            {
                I_q = I + mq*zs;
                gm[I_q] = bmdq[mq];
                gp[I_q] = -bpdq[mq]; 
                if (order[1] > 0)
                {
                    bpdq_trans[I_q] = bpdq[mq];
                    bmdq_trans[I_q] = bmdq[mq];                                                   
                }
            }

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

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

    maxcflblocks[blockIdx.z] = BlockReduce(temp_storage).Reduce(maxcfl,cub::Max());

    __syncthreads();


    /* ---------------------- Second order corrections and limiters --------------------*/  
    
    if (order[0] == 2)
    {
        ifaces_x = mx + 1;
        ifaces_y = my;
        num_ifaces = ifaces_x*ifaces_y;

        for(thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)

                        bpdq_trans[I_q] -= cqyy;      
                    } 
                }   
            }  
        }  
        __syncthreads();
