subroutine advecc2nd_corr_conservative(var, rhs)
! Removes the advection contribution from solid velocities, which should be
! close to zero but are not necessarily due to pressure correction.
! Has a fairly drastic effect on the initial flow, but the scalar is
! conserved throughout the simulation.
use modglobal, only : eps1, ib, ie, ih, jb, je, jh, kb, ke, kh, &
dx2i, dxi5, dy2i, dyi5, dzf, dzh2i, dzfi, dzhi, dzfi5
use modfields, only : u0, v0, w0
use modsubgriddata, only : ekh
use decomp_2d, only : zstart
real, intent(in) :: var(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh)
real, intent(inout) :: rhs(ib-ih:ie+ih,jb-jh:je+jh,kb :ke+kh)
integer :: i, j, k, n, m
do n = 1,bound_info_c%nbndptsrank
!n = bound_info_c%bndptsrank(m)
i = bound_info_c%bndpts_loc(n,1) - zstart(1) + 1
j = bound_info_c%bndpts_loc(n,2) - zstart(2) + 1
k = bound_info_c%bndpts_loc(n,3) - zstart(3) + 1
if ((abs(mask_u(i+1,j,k)) < eps1) .or. (abs(mask_c(i+1,j,k)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) + u0(i+1,j,k)*(var(i+1,j,k) + var(i,j,k))*dxi5
end if
if ((abs(mask_u(i,j,k)) < eps1) .or. (abs(mask_c(i-1,j,k)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) - u0(i,j,k)*(var(i-1,j,k) + var(i,j,k))*dxi5
end if
if ((abs(mask_v(i,j+1,k)) < eps1) .or. (abs(mask_c(i,j+1,k)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) + v0(i,j+1,k)*(var(i,j+1,k) + var(i,j,k))*dyi5
end if
if ((abs(mask_v(i,j,k)) < eps1) .or. (abs(mask_c(i,j-1,k)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) - v0(i,j,k)*(var(i,j-1,k) + var(i,j,k))*dyi5
end if
if ((abs(mask_w(i,j,k+1)) < eps1) .or. (abs(mask_c(i,j,k+1)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) + w0(i,j,k+1)*(var(i,j,k+1)*dzf(k) + var(i,j,k)*dzf(k+1))*dzhi(k+1)*dzfi5(k)
end if
if ((abs(mask_w(i,j,k)) < eps1) .or. (abs(mask_c(i,j,k-1)) < eps1)) then
rhs(i,j,k) = rhs(i,j,k) - w0(i,j,k)*(var(i,j,k-1)*dzf(k) + var(i,j,k)*dzf(k-1))*dzhi(k)*dzfi5(k)
end if
end do
end subroutine advecc2nd_corr_conservative