compute_terminal_health Subroutine

   subroutine compute_terminal_health(mesh, flow, params, fields, transport, energy, &
                                      projection_metric, kinetic_energy, tmin, tmax)
      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
      type(energy_fields_t), intent(in) :: energy
      real(rk), intent(out) :: projection_metric
      real(rk), intent(out) :: kinetic_energy
      real(rk), intent(out) :: tmin
      real(rk), intent(out) :: tmax

      integer :: c, lf, fid, ierr
      real(rk) :: vol, divu, flux_out, target, rho_cell
      real(rk) :: local_metric, local_ke, global_ke
      real(rk) :: local_temp_vec(2), global_temp_vec(2)
      logical :: have_temperature

      local_metric = zero
      local_ke = zero
      local_temp_vec = [huge(one), huge(one)]
      have_temperature = params%enable_energy .and. allocated(energy%T)

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

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

         divu = zero
         if (allocated(fields%face_flux)) then
            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
         end if

         target = zero
         if (params%enable_variable_density .and. allocated(fields%projection_divergence_source)) then
            if (c <= size(fields%projection_divergence_source)) target = fields%projection_divergence_source(c)
         end if
         local_metric = max(local_metric, abs(divu - target))

         rho_cell = params%rho
         if (params%enable_variable_density .and. allocated(transport%rho)) then
            if (c <= size(transport%rho)) rho_cell = max(transport%rho(c), tiny(one))
         end if
         local_ke = local_ke + 0.5_rk * rho_cell * vol * dot_product(fields%u(:, c), fields%u(:, c))

         if (have_temperature) then
            if (c <= size(energy%T)) then
               local_temp_vec(1) = min(local_temp_vec(1), energy%T(c))
               local_temp_vec(2) = min(local_temp_vec(2), -energy%T(c))
            end if
         end if
      end do

      call MPI_Allreduce(local_metric, projection_metric, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal projection health')

      call MPI_Allreduce(local_ke, global_ke, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal kinetic energy')
      kinetic_energy = global_ke

      if (have_temperature) then
         call MPI_Allreduce(local_temp_vec, global_temp_vec, 2, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing terminal temperature range')
         tmin = global_temp_vec(1)
         tmax = -global_temp_vec(2)
      else
         tmin = zero
         tmax = zero
      end if
   end subroutine compute_terminal_health