         {
             int I_q = I + mq*zs;  
             double gupdate = 0.5*dtdy*bpasdq[mq];
-            fm[I_q + 1] -= gupdate;
-            fp[I_q + 1] -= gupdate;
+            //fm[I_q + 1] -= gupdate;
+            //fp[I_q + 1] -= gupdate;
+            atomicAdd(&fm[I_q + 1],-gupdate);
+            atomicAdd(&fp[I_q + 1],-gupdate);
         }   
     } 
 
-    __syncthreads();
+    //__syncthreads();
 
     /* ------------------------------- Final update ----------------------------------- */
 

        {
            int I_q = I + mq*zs;  
            double gupdate = 0.5*dtdx*bmasdq[mq];
            gm[I_q - 1] -= gupdate;       
            gp[I_q - 1] -= gupdate;   
            //atomicAdd(&gm[I_q-1],-gupdate);
            //atomicAdd(&gp[I_q-1],-gupdate);
        }            
    }
    __syncthreads();


    /*     transverse-x

        {
            int I_q = I + mq*zs;  
            double gupdate = 0.5*dtdx*bpasdq[mq];
            //gm[I_q + ys] -= gupdate;
            //gp[I_q + ys] -= gupdate;
            atomicAdd(&gm[I_q + ys],-gupdate);
            atomicAdd(&gp[I_q + ys],-gupdate);
        }
        
    } 
    //__syncthreads();

    
    /* ----------------------------- Transverse : Y-faces ----------------------------- */
