do k=1,9
   e_k(1) = real(e((k-1)*2+1))
   e_k(2) = real(e((k-1)*2+2))
   DP = e_k(1)*U_x(i,j) + e_k(2)*U_y(i,j)
   f_eq = w(k)*rho(i,j)*(1.0+3.0*DP+(9.0/2.0)*(DP**2.0)-(3.0/2.0)*DPU)
   f(i,j,k)=f(i,j,k)-(1.0/tau)*(f(i,j,k)-f_eq)
enddo