function rho_np1 = rho_1st(vdx_fn, rho_n)
global delv; global delt;
[Nx,Nz] = size(rho_n);
rho_np1 = reshape(sum(-vdx_fn,2)*delv, Nx,Nz)*delt + rho_n;