subroutine diffc_corr(var, rhs, hi, hj, hk)
! Negate subgrid rhs contributions from solid points (added by diffc in modsubgrid)
use modglobal, only : eps1, ib, ie, jb, je, kb, ke, kh, &
dx2i, dxi5, dy2i, dyi5, dzf, dzh2i, dzfi, dzhi, dzfi5
use modsubgriddata, only : ekh
use decomp_2d, only : zstart
integer, intent(in) :: hi, hj, hk
real, intent(in) :: var(ib-hi:ie+hi,jb-hj:je+hj,kb-hk:ke+hk)
real, intent(inout) :: rhs(ib-hi:ie+hi,jb-hj:je+hj,kb :ke+hk)
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_c(i+1,j,k)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) - 0.5 * (ekh(i+1,j,k) + ekh(i,j,k)) * (var(i+1,j,k) - var(i,j,k))*dx2i
end if
if (abs(mask_c(i-1,j,k)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) + 0.5 * (ekh(i,j,k) + ekh(i-1,j,k)) * (var(i,j,k) - var(i-1,j,k))*dx2i
end if
if (abs(mask_c(i,j+1,k)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) - 0.5 * (ekh(i,j+1,k) + ekh(i,j,k)) * (var(i,j+1,k) - var(i,j,k))*dy2i
end if
if (abs(mask_c(i,j-1,k)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) + 0.5 * (ekh(i,j,k) + ekh(i,j-1,k)) * (var(i,j,k) - var(i,j-1,k))*dy2i
end if
if (abs(mask_c(i,j,k+1)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) - 0.5 * (dzf(k+1)*ekh(i,j,k) + dzf(k)*ekh(i,j,k+1)) &
* (var(i,j,k+1) - var(i,j,k))*dzh2i(k+1)*dzfi(k)
end if
if (abs(mask_c(i,j,k-1)) < eps1) then
rhs(i,j,k) = rhs(i,j,k) + 0.5 * (dzf(k-1)*ekh(i,j,k) + dzf(k)*ekh(i,j,k-1)) &
* (var(i,j,k) - var(i,j,k-1))*dzh2i(k)*dzfi(k)
end if
end do
end subroutine diffc_corr