             int y = subdomain.y0;
 
             int gi = IDX(x,y,z,lbm_sim->nx,lbm_sim->ny,lbm_sim->nz);
-            lbm_lattices_cuda_memcpy(h_lat, d_lat, gi, gi, ncc, cudaMemcpyDeviceToHost);
+            lbm_lattices_cuda_memcpy(lat0, lat1, gi, gi, ncc, kind);
         } else {
             // Read (#sny) partial planes
             for (int y = subdomain.y0; y <= subdomain.y1; y++) {
                 int gi = IDX(x,y,z,lbm_sim->nx,lbm_sim->ny,lbm_sim->nz);
-                lbm_lattices_cuda_memcpy(h_lat, d_lat, gi, gi, ncc, cudaMemcpyDeviceToHost);
+                lbm_lattices_cuda_memcpy(lat0, lat1, gi, gi, ncc, kind);
             }
         }
     } else {

    HANDLE_ERROR(cudaMemcpy(&lat0->bw[gi0], &lat1->bw[gi1], size, kind));
}

void lbm_lattices_read_subdomain(lbm_simulation* lbm_sim, lbm_lattices* h_lat, lbm_box_3d subdomain)
{
    lbm_lattices* d_lat = lbm_sim->switch_f0_f1 ? &lbm_sim->d_vars.f1 : &lbm_sim->d_vars.f0;
    size_t snx = std::abs(subdomain.x0 - subdomain.x1) + 1;
    size_t sny = std::abs(subdomain.y0 - subdomain.y1) + 1;
    size_t snz = std::abs(subdomain.z0 - subdomain.z1) + 1;

void lbm_lattices_write_subdomain(lbm_simulation* lbm_sim, lbm_lattices* h_lat, lbm_box_3d subdomain)
{
    lbm_lattices* d_lat = lbm_sim->switch_f0_f1 ? &lbm_sim->d_vars.f1 : &lbm_sim->d_vars.f0;
    lbm_lattices_cuda_memcpy_subdomain(lbm_sim, d_lat, h_lat, subdomain, cudaMemcpyHostToDevice);
}

void lbm_lattices_write(lbm_simulation* lbm_sim, lbm_lattices* h_lat)
