Doubts in Zou-He boundary condtion

Hello,
I am very new to the world of LBM. I am trying to simulate the benchmark problems with different boundary conditions.
I am confused with the algorithm of the Zou-He boundary condition. I want to use inlet velocity and zerogradient at outlet BC. Could you have a look at the following code and correct me where am I making mistake. Thank you

// 8 3 5 ^y
// |/ | x
// 2-0-1 —>
// /|
// 6 4 7

void LBM(int time) {

/// Swap populations
// The present code used old and new populations which are swapped at the beginning of each time step.
// This is sometimes called 'double-buffered' or 'ping-pong' algorithm.
// This way, the old populations are not overwritten during propagation.
// The resulting code is easier to write and to debug.
// The memory requirement for the populations is twice as large.

double ***swap_temp = pop_old;
pop_old = pop;
pop = swap_temp;

/// Lattice Boltzmann equation
// The lattice Boltzmann equation is solved in the following.
// The algorithm includes
// - computation of the lattice force
// - combined collision and propagation (faster than first collision and then propagation)

for (int X = 1; X < Nx-2; ++X) {
	for (int Y = 1; Y < Ny - 1; ++Y) {

		/// Compute equilibrium
		// The equilibrium populations are computed.
		
		equilibrium(density[X][Y], (velocity_x[X][Y] + (force_x[X][Y] + gravity) * tau / density[X][Y]), (velocity_y[X][Y] + (force_y[X][Y]) * tau / density[X][Y]));



		/// Compute new populations
		// This is the lattice Boltzmann equation (combined collision and propagation) including external forcing.
		

		pop[0][X][Y] = pop_old[0][X][Y] * (1 - omega) + pop_eq[0] * omega;
		pop[1][(X + 1)][Y] = pop_old[1][X][Y] * (1 - omega) + pop_eq[1] * omega;
		pop[2][X - 1][Y] = pop_old[2][X][Y] * (1 - omega) + pop_eq[2] * omega;
		pop[3][X][Y + 1] = pop_old[3][X][Y] * (1 - omega) + pop_eq[3] * omega;
		pop[4][X][Y - 1] = pop_old[4][X][Y] * (1 - omega) + pop_eq[4] * omega;
		pop[5][X + 1][Y + 1] = pop_old[5][X][Y] * (1 - omega) + pop_eq[5] * omega;
		pop[6][X - 1][Y - 1] = pop_old[6][X][Y] * (1 - omega) + pop_eq[6] * omega;
		pop[7][X + 1][Y - 1] = pop_old[7][X][Y] * (1 - omega) + pop_eq[7] * omega;
		pop[8][X - 1][Y + 1] = pop_old[8][X][Y] * (1 - omega) + pop_eq[8] * omega;


	}

//Bounce back condition for top and bottom wall
for (int X = 1; X < Nx-1; X++)
{

// Bottom wall (y = 0)


pop[3][X][1] = pop[4][X][0];
pop[5][X][1] = pop[6][X][0];
pop[8][X][1] = pop[7][X][0];

//Top wall (y = Ny - 1)

pop[4][X][Ny - 2] = pop[3][X][Ny - 1];
pop[6][X][Ny - 2] = pop[5][X][Ny - 1];
pop[7][X][Ny - 2] = pop[8][X][Ny - 1];

}
//Zou-He BC for inlet and outlet

for (int Y = 1; Y < Ny-1; Y++) {

//Left wall (X = 0)

pop[1][1][Y] = pop[2][1][Y] + (2 / 3)* inletdensity * inletvel;
pop[5][1][Y] = pop[6][1][Y] - 0.5*(pop[3][1][Y] - pop[4][1][Y]) + (1 / 6)*inletdensity* inletvel;
pop[7][1][Y] = pop[8][1][Y] + 0.5*(pop[3][1][Y] - pop[4][1][Y]) + (1 / 6)*inletdensity * inletvel;

//Right wall (X = Nx - 1)

pop[1][Nx - 1][Y] = pop[1][Nx - 2][Y];
pop[5][Nx - 1][Y] = pop[5][Nx - 2][Y];
pop[7][Nx - 1][Y] = pop[7][Nx - 2][Y];

}