energy_reconciliation_value Function

   real(rk) function energy_reconciliation_value(energy, transport, c, name) result(value)
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      integer, intent(in) :: c
      character(len=*), intent(in) :: name

      real(rk) :: rho_h_output, rho_h_operator, rho_h_recon, denom

      value = 0.0_rk
      rho_h_output = 0.0_rk
      rho_h_operator = 0.0_rk
      rho_h_recon = 0.0_rk

      if (.not. allocated(energy%h)) return
      if (.not. allocated(transport%rho)) return
      if (c < 1 .or. c > size(energy%h) .or. c > size(transport%rho)) return

      rho_h_output = transport%rho(c) * energy%h(c)

      ! Before the first energy update, no operator-consistent cellwise state
      ! exists.  Fall back to the output-state value so step-0 visualization has
      ! a defined, zero reconciliation field.
      rho_h_operator = rho_h_output
      if (energy%operator_consistent_rho_h_available == 1 .and. &
          allocated(energy%operator_consistent_rho_h)) then
         if (size(energy%operator_consistent_rho_h) >= c) then
            rho_h_operator = energy%operator_consistent_rho_h(c)
         end if
      end if

      rho_h_recon = rho_h_output - rho_h_operator

      select case (trim(name))
      case ('rho_h_output_state')
         value = rho_h_output
      case ('rho_h_operator_consistent')
         value = rho_h_operator
      case ('rho_h_density_reconciliation')
         value = rho_h_recon
      case ('relative_rho_h_density_reconciliation')
         denom = max(abs(rho_h_output), abs(rho_h_operator), tiny(1.0_rk))
         value = abs(rho_h_recon) / denom
      case default
         value = 0.0_rk
      end select

   end function energy_reconciliation_value