   gx = 0.0;
   gy = 0.0;
   gz = 0.0;
+  
+  #ifdef GRAVITY_5_POINTS_GRADIENT
+  int id_ll, id_rr;
+  Real pot_ll, pot_rr;
+  #endif
   
   #ifdef COUPLE_DELTA_E_KINETIC
   Real Ekin_0, Ekin_1;

    id_r = (xid+1) + (yid)*nx + (zid)*nx*ny;
    pot_l = dev_potential[id_l];
    pot_r = dev_potential[id_r];
    gx = -0.5*( pot_r - pot_l ) / dx;
    
    //Get Y componet of gravity field
    id_l = (xid) + (yid-1)*nx + (zid)*nx*ny;
    id_r = (xid) + (yid+1)*nx + (zid)*nx*ny;
    pot_l = dev_potential[id_l];
    pot_r = dev_potential[id_r];
    gy = -0.5*( pot_r - pot_l ) / dy;
    
    //Get Z componet of gravity field
    id_l = (xid) + (yid)*nx + (zid-1)*nx*ny;
    id_r = (xid) + (yid)*nx + (zid+1)*nx*ny;
    pot_l = dev_potential[id_l];
    pot_r = dev_potential[id_r];
    gz = -0.5*( pot_r - pot_l ) / dz;
    
    
    dev_conserved[  n_cells + id] += 0.5*dt*gx*(d + d_n);

    if ( dev_conserved[            id] < density_floor ){
      if (dev_conserved[            id] > 0){  
        dens_0 = dev_conserved[            id];
        // printf("###Thread density change  %f -> %f \n", dens_0, density_floor );
        dev_conserved[            id] = density_floor;
        // Scale the conserved values to the new density
        dev_conserved[1*n_cells + id] *= (density_floor / dens_0);
