compute_projection_audit_cell_balance Subroutine

   subroutine compute_projection_audit_cell_balance(mesh, fields, c, div_c, mdot_div_c, &
                                                    boundary_face_count, physical_boundary_face_count)
      type(mesh_t), intent(in) :: mesh
      type(flow_fields_t), intent(in) :: fields
      integer, intent(in) :: c
      real(rk), intent(out) :: div_c
      real(rk), intent(out) :: mdot_div_c
      integer, intent(out) :: boundary_face_count
      integer, intent(out) :: physical_boundary_face_count

      integer :: lf, fid, nb
      real(rk) :: vol, flux_out, mdot_out, mass_flux_value

      div_c = 0.0_rk
      mdot_div_c = 0.0_rk
      boundary_face_count = 0
      physical_boundary_face_count = 0

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

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

         nb = mesh%faces(fid)%neighbor
         if (nb <= 0) then
            boundary_face_count = boundary_face_count + 1
            if (mesh%faces(fid)%periodic_neighbor <= 0) physical_boundary_face_count = physical_boundary_face_count + 1
         end if

         if (allocated(fields%face_flux)) then
            if (mesh%faces(fid)%owner == c) then
               flux_out = fields%face_flux(fid)
            else
               flux_out = -fields%face_flux(fid)
            end if
            div_c = div_c + flux_out / vol
         end if

         if (allocated(fields%mass_flux)) then
            mass_flux_value = fields%mass_flux(fid)
            if (mesh%faces(fid)%owner == c) then
               mdot_out = mass_flux_value
            else
               mdot_out = -mass_flux_value
            end if
            mdot_div_c = mdot_div_c + mdot_out / vol
         end if
      end do

   end subroutine compute_projection_audit_cell_balance