         {
             int I_q = I + mq*zs;  
             double gupdate = 0.5*dtdy*bmasdq[mq];
-            //fm[I_q] -= gupdate;        
-            //fp[I_q] -= gupdate;
+#if CUDA_COMPUTE==3            
+            fm[I_q] -= gupdate;        
+            fp[I_q] -= gupdate;
+#else            
             atomicAdd(&fm[I_q], -gupdate);        
             atomicAdd(&fp[I_q], -gupdate);
 
             gupdate = 0.5*dtdy*bpasdq[mq];
             atomicAdd(&fm[I_q + 1],-gupdate);
             atomicAdd(&fp[I_q + 1],-gupdate);
+#endif            
 
         }
     }

        {
            int I_q = I + mq*zs;  
            double gupdate = 0.5*dtdy*bmasdq[mq];
            //fm[I_q - ys] -= gupdate;        
            //fp[I_q - ys] -= gupdate;
            atomicAdd(&fm[I_q - ys],-gupdate);        
            atomicAdd(&fp[I_q - ys],-gupdate);

            gupdate = 0.5*dtdy*bpasdq[mq];
            atomicAdd(&fm[I_q - ys + 1],-gupdate);
            atomicAdd(&fp[I_q - ys + 1],-gupdate);                
        }
    }
    __syncthreads();

             |     |     
    */              

#if CUDA_COMPUTE==3
    for(int thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
    { 
        int ix = thread_index % ifaces_x;
