compute_global_observables Subroutine

   subroutine compute_global_observables(mesh, flow, fields, species, transport, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(in) :: species
      type(transport_properties_t), intent(in) :: transport
      type(solver_stats_t), intent(inout) :: stats

      integer :: c
      real(rk) :: local_max_vel, local_min_y, local_total_mass
      real(rk) :: global_max_vel, global_min_y, global_total_mass
      real(rk) :: u, v, w, vel_mag

      local_max_vel = zero
      local_min_y = huge(one)
      local_total_mass = zero

      do c = 1, mesh%ncells
         if (.not. flow%owned(c)) cycle

         u = fields%u(1, c)
         v = fields%u(2, c)
         w = fields%u(3, c)
         vel_mag = sqrt(u**2 + v**2 + w**2)
         if (vel_mag > local_max_vel) local_max_vel = vel_mag

         local_total_mass = local_total_mass + (transport%rho(c) * mesh%cells(c)%volume)

         if (species%nspecies > 0) then
            local_min_y = min(local_min_y, minval(species%Y(:, c)))
         end if
      end do

      call MPI_Allreduce(local_max_vel, global_max_vel, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm)
      call MPI_Allreduce(local_total_mass, global_total_mass, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm)
      
      if (species%nspecies > 0) then
         call MPI_Allreduce(local_min_y, global_min_y, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm)
      else
         global_min_y = zero
      end if

      stats%max_velocity = global_max_vel
      stats%total_mass = global_total_mass
      stats%min_species_y = global_min_y

   end subroutine compute_global_observables