             // boundaries
             cudaMemcpy(Un_h, Un_d, nx*ny*nlayers*4*sizeof(float), cudaMemcpyDeviceToHost);
             //bcs_fv(Un_h, nx, ny, nlayers, ng);
-            bcs_mpi(Un_h, nx, ny, nlayers, ng, comm, status, rank, n_processes);
+            if (n_processes == 1) {
+                bcs_fv(Un_h, nx, ny, nlayers, ng);
+            } else {
+                bcs_mpi(Un_h, nx, ny, nlayers, ng, comm, status, rank, n_processes);
+            }
             cudaMemcpy(Un_d, Un_h, nx*ny*nlayers*4*sizeof(float), cudaMemcpyHostToDevice);
 
             cudaError_t err = cudaGetLastError();

    // copy back flux
    cudaMemcpy(F_h, F_d, nx*ny*nlayers*4*sizeof(float), cudaMemcpyDeviceToHost);
    //bcs_fv(F_h, nx, ny, nlayers, ng);
    bcs_mpi(F_h, nx, ny, nlayers, ng, comm, status, rank, n_processes);

    for (int y = 0; y < ny; y++) {
        for (int x = 0; x < nx; x++) {


        // easiest (but maybe not most balanced way) is to split kernels into strips
        // not enough kernels to fill all processes. This would be inefficient if ny > nx, so need to fix this to look in the y direction if this is the case.
        if (kernels_y < n_processes) {
            for (int i = 0; i < kernels_y; i++) {
                kernels[i].x = kernels_x;
                kernels[i].y = 1;
            }
            // blank out the other ones
            for (int i = kernels_y; i < n_processes; i++) {
                kernels[i].x = 0;
                kernels[i].y = 0;
            }
