subroutine fluxtop(field, ek, flux)
use modglobal, only:ib, ie, ih, jb, je, jh, kb, ke, kh, dzf, dzh, dzhi, eps1
real, intent(inout) :: field(ib - ih:ie + ih, jb - jh:je + jh, kb - kh:ke + kh)
real, intent(in) :: ek(ib - ih:ie + ih, jb - jh:je + jh, kb - kh:ke + kh)
real, intent(in) :: flux
!
if (abs(flux) .le. eps1) then !it's zero-flux, we don't need to do the calculation
field(:, :, ke + 1) = field(:, :, ke)
else
field(:, :, ke + 1) = field(:, :, ke) + dzh(ke + 1)*flux/(dzhi(ke + 1)*(0.5*(dzf(ke)*ek(:, :, ke + 1) + dzf(ke + 1)*ek(:, :, ke))))
end if
!
end subroutine fluxtop