lowmach_cell_u_dot_grad_rho Function

   real(rk) function lowmach_cell_u_dot_grad_rho(mesh, fields, transport, c) result(value)
      type(mesh_t), intent(in) :: mesh
      type(flow_fields_t), intent(in) :: fields
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: c

      integer :: lf, fid, nb
      real(rk) :: vol, flux_out, rho_c, rho_nb, rho_f

      value = 0.0_rk
      if (.not. allocated(fields%face_flux)) return
      if (.not. allocated(transport%rho)) return
      if (c < 1 .or. c > mesh%ncells) return
      if (size(transport%rho) < c) return

      vol = mesh%cells(c)%volume
      if (vol <= 0.0_rk) return

      rho_c = transport%rho(c)

      do lf = 1, mesh%ncell_faces(c)
         fid = mesh%cell_faces(lf, c)

         if (mesh%faces(fid)%owner == c) then
            flux_out = fields%face_flux(fid)
            nb = mesh%faces(fid)%neighbor
            if (nb <= 0) nb = mesh%faces(fid)%periodic_neighbor
         else
            flux_out = -fields%face_flux(fid)
            nb = mesh%faces(fid)%owner
         end if

         if (nb > 0 .and. nb <= mesh%ncells .and. size(transport%rho) >= nb) then
            rho_nb = transport%rho(nb)
            rho_f = 0.5_rk * (rho_c + rho_nb)
         else
            rho_f = rho_c
         end if

         value = value + flux_out * (rho_f - rho_c)
      end do

      value = value / vol

   end function lowmach_cell_u_dot_grad_rho