         double *const qr     = start;          /* meqn   */
         double *const qd     = qr + meqn;      /* meqn   */
         double *const bmdq   = qd + meqn;    /* meqn   */
-        double *const bpdq   = bmdq + meqn;    /* meqn   */
+        double *const aux1   = bmdq + meqn;    /* 2*maux */
+        double *const aux2   = aux1 + 2*maux;  /* 2*maux */
+        double *const aux3   = aux2 + 2*maux;  /* 2*maux */
+        double *const bmasdq = aux3 + 2*maux;  /* meqn   */
+        double *const bpasdq = bmasdq + meqn;  /* meqn   */
+        for(int mq = 0; mq < meqn; mq++)
+        {
+            int I_q = I + mq*zs;
+
+            qr[mq] = qold[I_q];
+            qd[mq] = qold[I_q - ys];
+
+            bmdq[mq] = bmdq_trans[I_q];
+            //bpdq[mq] = bpdq_trans[I_q];
+        }            
+
+        for(int imp = 0; imp < 2; imp++)
+        {
+            for(int m = 0; m < maux; m++)
+            {
+                int I_aux = I + m*zs;
+                int k = imp*maux + m;
+                aux1[k] = aux[I_aux - 1 + ys*(imp - 1)];
+                aux2[k] = aux[I_aux     + ys*(imp - 1)];
+                aux3[k] = aux[I_aux + 1 + ys*(imp - 1)];
+            }
+        }
+
+        rpt2(1,meqn,mwaves,maux,qd,qr,aux1,aux2,aux3,0,bmdq,bmasdq,bpasdq);
+
+        for(int mq = 0; mq < meqn; mq++)
+        {
+            int I_q = I + mq*zs;  
+            double gupdate = 0.5*dtdy*bmasdq[mq];
+            fm[I_q - ys] -= gupdate;        
+            fp[I_q - ys] -= gupdate;
+        }
+
+    }
+    __syncthreads();
+
+
+    /*  transverse-y
+    
+             |     |     
+        -----|-----|-----
+             |     |     
+             |  q  |      
+             |     |     
+        -----|-XXX-|-----
+             |   v |     
+             0   --0     
+             |     |     
+        -----|-----|-----
+             |     |     
+    */              
+
+    for(int thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
+    { 
+        int ix = thread_index % ifaces_x;
+        int iy = thread_index/ifaces_x;
+
+        int I =  (iy + mbc)*ys + (ix + mbc-1);
+
+
+        double *const qr     = start;          /* meqn   */
+        double *const qd     = qr + meqn;      /* meqn   */
+        double *const bmdq   = qd + meqn;    /* meqn   */
+        double *const aux1   = bmdq + meqn;    /* 2*maux */
+        double *const aux2   = aux1 + 2*maux;  /* 2*maux */
+        double *const aux3   = aux2 + 2*maux;  /* 2*maux */
+        double *const bmasdq = aux3 + 2*maux;  /* meqn   */
+        double *const bpasdq = bmasdq + meqn;  /* meqn   */
+        for(int mq = 0; mq < meqn; mq++)
+        {
+            int I_q = I + mq*zs;
+
+            qr[mq] = qold[I_q];
+            qd[mq] = qold[I_q - ys];
+
+            bmdq[mq] = bmdq_trans[I_q];
+            //bpdq[mq] = bpdq_trans[I_q];
+        }            
+
+        for(int imp = 0; imp < 2; imp++)
+        {
+            for(int m = 0; m < maux; m++)
+            {
+                int I_aux = I + m*zs;
+                int k = imp*maux + m;
+                aux1[k] = aux[I_aux - 1 + ys*(imp - 1)];
+                aux2[k] = aux[I_aux     + ys*(imp - 1)];
+                aux3[k] = aux[I_aux + 1 + ys*(imp - 1)];
+            }
+        }
+
+        rpt2(1,meqn,mwaves,maux,qd,qr,aux1,aux2,aux3,0,bmdq,bmasdq,bpasdq);
+
+        for(int mq = 0; mq < meqn; mq++)
+        {
+            int I_q = I + mq*zs;  
+            int gupdate = 0.5*dtdy*bpasdq[mq];
+            fm[I_q - ys + 1] -= gupdate;
+            fp[I_q - ys + 1] -= gupdate;                
+        }
+    }
+    __syncthreads();
+
+
+
+    /*  transverse-y
+    
+             |     |     
+        -----|-----|-----
+             |  q  |     
+             O--   |           
+             | ^   |     
+        -----|-XXX-|-----
+             |     |     
+             |     | 
+             |     |     
+        -----|-----|-----
+             |     |     
+    */              
+
+    for(int thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
+    { 
+        int ix = thread_index % ifaces_x;
+        int iy = thread_index/ifaces_x;
+
+        int I =  (iy + mbc)*ys + (ix + mbc-1);
+
+
+        double *const qr     = start;          /* meqn   */
+        double *const qd     = qr + meqn;      /* meqn   */
+        double *const bpdq   = qd + meqn;    /* meqn   */
         double *const aux1   = bpdq + meqn;    /* 2*maux */
         double *const aux2   = aux1 + 2*maux;  /* 2*maux */
         double *const aux3   = aux2 + 2*maux;  /* 2*maux */

                int I_aux = I + m*zs;
                auxr[m] = aux[I_aux];
                auxl[m] = aux[I_aux - 1];
                //auxd[m] = aux[I_aux - ys];
            }               

            

            ql[mq] = qold[I_q-1];
            qr[mq] = qold[I_q];

            //amdq[mq] = amdq_trans[I_q];
            apdq[mq] = apdq_trans[I_q];
        }            

