                 double fout_ne, fout_e, fout_se, fout_n, fout_c, fout_s, fout_nw, fout_w, fout_sw,
                        fout_te, fout_tn, fout_tc, fout_ts, fout_tw, 
                        fout_be, fout_bn, fout_bc, fout_bs, fout_bw;
+                double feq_ne, feq_e, feq_se, feq_n, feq_c, feq_s, feq_nw, feq_w, feq_sw, 
+                       feq_te, feq_tn, feq_tc, feq_ts, feq_tw, 
+                       feq_be, feq_bn, feq_bc, feq_bs, feq_bw;
 
                 fin_ne = f0.ne[gi];
                 fin_e  = f0.e [gi];

    lbm_sim->dimComputationGrid = dimComputationGrid;
    lbm_sim->dimComputationBlock = dimComputationBlock;

    dim3 dimRightWallGrid(1, min((unsigned long)ny, (unsigned long)prop.maxGridSize[1]));
    dim3 dimRightWallBlock(1);
    lbm_sim->dimRightWallGrid = dimRightWallGrid;
    lbm_sim->dimRightWallBlock = dimRightWallBlock;

    lbm_sim->shared_mem_size = 6 * sizeof(double) * BLOCK_SIZE;


    return lbm_sim;
}

#define OMEGA    ((double)1. / (3*(NULB)+0.5))  // Relaxation parameter

#define BLOCK_SIZE 64

#ifndef SHARED_MEM_SIZE
#define SHARED_MEM_SIZE 0
#endif

#define SQUARE(a) ((a)*(a))
#define GPU_SQUARE(a) __dmul_rn(a,a)
