write_variable_density_energy_debug Subroutine

   subroutine write_variable_density_energy_debug(mesh, flow, params, fields, energy, label, species_Y, transport, step)
      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(energy_fields_t), intent(in) :: energy
      character(len=*), intent(in) :: label
      real(rk), intent(in), optional :: species_Y(:,:)
      type(transport_properties_t), intent(in), optional :: transport
      integer, intent(in), optional :: step

      integer :: c, k, ierr
      integer :: local_bad_cell
      real(rk) :: local_min_h, local_max_h, local_min_T, local_max_T
      real(rk) :: local_min_rho, local_max_rho, local_min_rho_old, local_max_rho_old
      real(rk) :: local_min_S, local_max_S
      real(rk) :: global_min_h, global_max_h, global_min_T, global_max_T
      real(rk) :: global_min_rho, global_max_rho, global_min_rho_old, global_max_rho_old
      real(rk) :: global_min_S, global_max_S
      real(rk) :: local_min_vec(5), local_max_vec(5)
      real(rk) :: global_min_vec(5), global_max_vec(5)
      logical :: have_transport, suspicious

      if (.not. params%enable_variable_density) return
      if (.not. params%variable_density_debug) return
      if (present(step)) then
         if (mod(step, max(1, params%variable_density_debug_interval)) /= 0) return
      end if
      if (.not. allocated(energy%h) .or. .not. allocated(energy%T)) return

      associate(dummy_mesh_ncells => mesh%ncells)
      end associate

      have_transport = present(transport)

      local_min_h = huge(zero)
      local_max_h = -huge(zero)
      local_min_T = huge(zero)
      local_max_T = -huge(zero)
      local_min_rho = huge(zero)
      local_max_rho = -huge(zero)
      local_min_rho_old = huge(zero)
      local_max_rho_old = -huge(zero)
      local_min_S = huge(zero)
      local_max_S = -huge(zero)
      local_bad_cell = -1

      do c = flow%first_cell, flow%last_cell
         local_min_h = min(local_min_h, energy%h(c))
         local_max_h = max(local_max_h, energy%h(c))
         local_min_T = min(local_min_T, energy%T(c))
         local_max_T = max(local_max_T, energy%T(c))

         if (have_transport) then
            if (allocated(transport%rho)) then
               local_min_rho = min(local_min_rho, transport%rho(c))
               local_max_rho = max(local_max_rho, transport%rho(c))
            end if
            if (allocated(transport%rho_old)) then
               local_min_rho_old = min(local_min_rho_old, transport%rho_old(c))
               local_max_rho_old = max(local_max_rho_old, transport%rho_old(c))
            end if
         end if

         if (allocated(fields%divergence_source)) then
            local_min_S = min(local_min_S, fields%divergence_source(c))
            local_max_S = max(local_max_S, fields%divergence_source(c))
         end if

         if (local_bad_cell < 0) then
            if (energy%h(c) < -1.0e5_rk .or. energy%h(c) > 1.0e7_rk) then
               local_bad_cell = c
            end if
         end if
      end do

      if (.not. have_transport) then
         local_min_rho = zero
         local_max_rho = zero
         local_min_rho_old = zero
         local_max_rho_old = zero
      end if
      if (.not. allocated(fields%divergence_source)) then
         local_min_S = zero
         local_max_S = zero
      end if

      local_min_vec = [local_min_h, local_min_T, local_min_rho, local_min_rho_old, local_min_S]
      local_max_vec = [local_max_h, local_max_T, local_max_rho, local_max_rho_old, local_max_S]

      call MPI_Allreduce(local_min_vec, global_min_vec, 5, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing vd debug minima')
      call MPI_Allreduce(local_max_vec, global_max_vec, 5, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing vd debug maxima')

      global_min_h = global_min_vec(1)
      global_min_T = global_min_vec(2)
      global_min_rho = global_min_vec(3)
      global_min_rho_old = global_min_vec(4)
      global_min_S = global_min_vec(5)
      global_max_h = global_max_vec(1)
      global_max_T = global_max_vec(2)
      global_max_rho = global_max_vec(3)
      global_max_rho_old = global_max_vec(4)
      global_max_S = global_max_vec(5)

      if (flow%rank == 0) then
         write(output_unit,'(a,a,a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5),a,2(1x,es13.5))') &
            'vd-energy-debug[', trim(label), '] h=', global_min_h, global_max_h, &
            ' T=', global_min_T, global_max_T, &
            ' rho=', global_min_rho, global_max_rho, &
            ' rho_old=', global_min_rho_old, global_max_rho_old, &
            ' S=', global_min_S, global_max_S
      end if

      suspicious = local_bad_cell > 0
      if (suspicious) then
         c = local_bad_cell
         write(output_unit,'(a,i0,a,i0,a,a,a,es16.8,a,es16.8)') &
            'vd-energy-bad-cell rank=', flow%rank, ' cell=', c, ' label=', trim(label), &
            ' h=', energy%h(c), ' T=', energy%T(c)

         if (have_transport) then
            if (allocated(transport%rho) .and. allocated(transport%rho_old)) then
               write(output_unit,'(a,es16.8,a,es16.8)') &
                  '  rho=', transport%rho(c), ' rho_old=', transport%rho_old(c)
            end if
         end if

         if (allocated(fields%divergence_source)) then
            write(output_unit,'(a,es16.8)') '  divergence_source=', fields%divergence_source(c)
         end if

         if (present(species_Y)) then
            do k = 1, min(params%nspecies, size(species_Y, 1))
               write(output_unit,'(a,i0,a,es16.8)') '  Y(', k, ')=', species_Y(k, c)
            end do
         end if
      end if

   end subroutine write_variable_density_energy_debug