             /* Start at first non-ghost interior cell */
             int I = (iy + mbc)*ys + ix + mbc - 1;
 
+            double *s = start;
+            double *wave = s + mwaves;
             for(int mw = 0; mw < mwaves; mw++)
             {
                 /* ------------------------------- Y-directions --------------------------- */

    
    */              


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


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

        {
            /* ------------------------ Normal solve in X direction ------------------- */
            double *const ql     = start;                 /* meqn        */
            double *const qr     = ql     + meqn;         /* meqn        */
            double *const auxl   = qr     + meqn;         /* maux        */
            double *const auxr   = auxl   + maux;         /* maux        */
            double *const s      = auxr   + maux;         /* mwaves      */
            double *const wave   = s      + mwaves;       /* meqn*mwaves */
            double *const amdq   = wave   + meqn*mwaves;  /* meqn        */
            double *const apdq   = amdq   + meqn;         /* meqn        */
            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];
            }               

            
            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 ------------------- */
            double *const ql     = start;                 /* meqn        */
            double *const qd     = ql     + meqn;         /* meqn        */
            double *const auxl   = qd     + meqn;         /* maux        */
            double *const auxd   = auxl   + maux;         /* maux        */
            double *const s      = auxr   + maux;         /* mwaves      */
            double *const wave   = s      + mwaves;       /* meqn*mwaves */
            double *const bmdq   = wave   + meqn*mwaves;  /* meqn        */
            double *const bpdq   = amdq   + meqn;         /* meqn        */

            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];
                    }
                }
            }
        }
