         }      
         __syncthreads(); /* Needed to be sure all aux variables are available below */ 
     } 
+#endif    
 
     /* -------------------------- Compute fluctuations -------------------------------- */
 


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

        for(int mq = 0; mq < meqn; mq++)
        {
            int 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(int m = 0; m < maux; m++)
        {
            int I_aux = I + m*zs;
            auxr[m] = aux[I_aux];
            auxl[m] = aux[I_aux - 1];
            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 (int mq = 0; mq < meqn; mq++) 
        {
            int 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(int mw = 0; mw < mwaves; mw++)
        {
            maxcfl = max(maxcfl,abs(s[mw]*dtdx));

            if (order[0] == 2)
            {                    
                int I_speeds = I + mw*zs;
                speeds[I_speeds] = s[mw];
                for(int mq = 0; mq < meqn; mq++)
                {
                    int k = mw*meqn + mq;
                    int 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 (int mq = 0; mq < meqn; mq++) 
        {
            int I_q = I + mq*zs;
            gm[I_q] = bmdq[mq];
            gp[I_q] = -bpdq[mq]; 
            if (order[1] > 0)
            {
                bmdq_trans[I_q] = bmdq[mq];                                                   
                bpdq_trans[I_q] = bpdq[mq];
            }
        }

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

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


    double* start  = shared_mem + mwork*threadIdx.x;

#if 0
    double *const ql     = start;                 /* meqn        */
    double *const qr     = ql     + meqn;         /* meqn        */
    double *const qd     = qr     + meqn;         /* meqn        */
