compute_variable_density_console_diagnostics Subroutine

   subroutine compute_variable_density_console_diagnostics(mesh, flow, params, fields, transport, time, &
                                                           max_divu, max_s_projection, &
                                                           max_lm_res_projection, max_lm_res_current, &
                                                           mass_balance_defect)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(transport_properties_t), intent(in) :: transport
      real(rk), intent(in) :: time
      real(rk), intent(out) :: max_divu
      real(rk), intent(out) :: max_s_projection
      real(rk), intent(out) :: max_lm_res_projection
      real(rk), intent(out) :: max_lm_res_current
      real(rk), intent(out) :: mass_balance_defect

      integer :: c, lf, fid, f, owner, nb, ierr
      real(rk) :: vol, divu, flux_out, s_current, s_projection
      real(rk) :: local_max_divu, local_max_s_projection
      real(rk) :: local_max_lm_res_projection, local_max_lm_res_current
      real(rk) :: local_mass, global_mass, local_net_mass_flux, net_mass_flux
      real(rk) :: dt_since_previous, mass_rate
      real(rk) :: local_max_vec(4), global_max_vec(4)
      real(rk) :: local_sum_vec(2), global_sum_vec(2)
      logical, save :: initialized = .false.
      real(rk), save :: previous_mass = 0.0_rk
      real(rk), save :: previous_time = 0.0_rk

      max_divu = 0.0_rk
      max_s_projection = 0.0_rk
      max_lm_res_projection = 0.0_rk
      max_lm_res_current = 0.0_rk
      mass_balance_defect = 0.0_rk

      if (.not. params%enable_variable_density) return
      if (.not. allocated(fields%face_flux)) return

      local_max_divu = 0.0_rk
      local_max_s_projection = 0.0_rk
      local_max_lm_res_projection = 0.0_rk
      local_max_lm_res_current = 0.0_rk
      local_mass = 0.0_rk
      local_net_mass_flux = 0.0_rk

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

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

         divu = 0.0_rk
         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)
            else
               flux_out = -fields%face_flux(fid)
            end if
            divu = divu + flux_out / vol
         end do

         s_current = 0.0_rk
         if (allocated(fields%divergence_source)) then
            if (size(fields%divergence_source) >= c) s_current = fields%divergence_source(c)
         end if

         s_projection = s_current
         if (allocated(fields%projection_divergence_source)) then
            if (size(fields%projection_divergence_source) >= c) s_projection = fields%projection_divergence_source(c)
         end if

         local_max_divu = max(local_max_divu, abs(divu))
         local_max_s_projection = max(local_max_s_projection, abs(s_projection))
         local_max_lm_res_projection = max(local_max_lm_res_projection, abs(divu - s_projection))
         local_max_lm_res_current = max(local_max_lm_res_current, abs(divu - s_current))

         if (allocated(transport%rho)) then
            if (size(transport%rho) >= c) local_mass = local_mass + transport%rho(c) * vol
         end if
      end do

      if (allocated(fields%mass_flux)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0) then
               if (owner >= 1 .and. owner <= mesh%ncells) then
                  if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + fields%mass_flux(f)
               end if
            end if
         end do
      else if (allocated(fields%face_flux) .and. allocated(transport%rho)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0) then
               if (owner >= 1 .and. owner <= mesh%ncells) then
                  if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + transport%rho(owner) * fields%face_flux(f)
               end if
            end if
         end do
      end if

      local_max_vec = [local_max_divu, local_max_s_projection, &
                       local_max_lm_res_projection, local_max_lm_res_current]
      call MPI_Allreduce(local_max_vec, global_max_vec, 4, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density max diagnostics')
      max_divu = global_max_vec(1)
      max_s_projection = global_max_vec(2)
      max_lm_res_projection = global_max_vec(3)
      max_lm_res_current = global_max_vec(4)

      local_sum_vec = [local_mass, local_net_mass_flux]
      call MPI_Allreduce(local_sum_vec, global_sum_vec, 2, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density sum diagnostics')
      global_mass = global_sum_vec(1)
      net_mass_flux = global_sum_vec(2)

      if (initialized) then
         dt_since_previous = time - previous_time
         if (dt_since_previous > tiny(1.0_rk)) then
            mass_rate = (global_mass - previous_mass) / dt_since_previous
            mass_balance_defect = mass_rate + net_mass_flux
         else
            mass_balance_defect = 0.0_rk
         end if
      else
         mass_balance_defect = 0.0_rk
         initialized = .true.
      end if

      previous_mass = global_mass
      previous_time = time

   end subroutine compute_variable_density_console_diagnostics