     for (int j = 0; j < kernels[rank].y; j++) {
        kx_offset = 0;
        for (int i = 0; i < kernels[rank].x; i++) {
-            //cout << "compressible_from_swe\n";
             swe_from_compressible<<<blocks[k_offset + j * kernels[rank].x + i], threads[k_offset + j * kernels[rank].x + i]>>>(q_fd, qf_swe, nxf, nyf, gamma_up, rho, gamma, kx_offset, ky_offset);
             kx_offset += blocks[k_offset + j * kernels[rank].x + i].x *
                 threads[k_offset + j * kernels[rank].x + i].x - 2*ng;


        // fluxes

        //printf("x, y: %i, %i, q_p: %f, Un_d: %f\n", x, y, q_p[0], Un_d[offset+0]);

        flux_func(q_p, f, true, gamma_up_d, alpha, beta_d, gamma);

        //printf("q_p, f: %f, %f\n", q_p[0], f[0]);

        //printf("Called flux func?\n");

        for (int i = 0; i < vec_dim; i++) {
            qx_plus_half[offset + i] = q_p[i];
            fx_plus_half[offset + i] = f[i];
        }

        //printf("fx_p_d (%i, %i): %f\n", x, y, fx_plus_half[offset]);

        flux_func(q_m, f, true, gamma_up_d, alpha, beta_d, gamma);


    int y = ky_offset + blockIdx.y * blockDim.y + threadIdx.y;
    int offset = y * nx + x;

    if ((x < nx) && (y < ny)) {
        float * q_swe;
        q_swe = (float *)malloc(3 * sizeof(float));

        for (int i = 0; i < 3; i++) {
            q_swe[i] = q[offset * 3 + i];
        }

