         ix = thread_index % ifaces_x;
         iy = thread_index/ifaces_y;
 
-        I = (iy + mbc + 1)*ys+ (ix + mbc);
+        I = (iy + mbc)*ys+ (ix + mbc + 1);
 
         //if (0 < ix && ix < mx + 1 && iy < my)   /* Is this needed? */
         {

            ix = thread_index % ifaces_x;
            iy = thread_index/ifaces_y;

            /* Start at first interior cell */
            I = (ix + mbc)*xs + (iy + mbc)*ys;

            //if (ix < mx + 1 && iy < my + 1)   /* Is this needed? */

        __syncthreads(); /* Needed to be sure all aux variables are available below */ 
    } 

    /* -------------------------- Compute fluctuations -------------------------------- */

    for(thread_index = threadIdx.x; thread_index < num_ifaces; thread_index += blockDim.x)
    {
        ix = thread_index % ifaces_x;
        iy = thread_index/ifaces_y;

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

        //if (ix < mx + 2*mbc-1 && iy < my + 2*mbc-1)
        {
