allocate_energy Subroutine

   subroutine allocate_energy(mesh, energy)
      type(mesh_t), intent(in) :: mesh
      type(energy_fields_t), intent(inout) :: energy

      call finalize_energy(energy)

      allocate(energy%T(mesh%ncells))
      allocate(energy%T_old(mesh%ncells))
      allocate(energy%h(mesh%ncells))
      allocate(energy%h_old(mesh%ncells))
      allocate(energy%rhs_old(mesh%ncells))
      allocate(energy%qrad(mesh%ncells))
      allocate(energy%cp(mesh%ncells))
      allocate(energy%lambda(mesh%ncells))
      allocate(energy%rho_thermo(mesh%ncells))
      allocate(energy%operator_consistent_rho_h(mesh%ncells))
      allocate(energy%species_enthalpy_diffusion(mesh%ncells))
      allocate(energy%chem_heat_release_rate(mesh%ncells))
      allocate(energy%chem_dT_dt(mesh%ncells))
      allocate(energy%chem_dY_max_dt(mesh%ncells))
      allocate(energy%chem_rel_rho_change_dt(mesh%ncells))
      allocate(energy%chem_source_alpha(mesh%ncells))
      allocate(energy%chem_active(mesh%ncells))

      energy%T = zero
      energy%T_old = zero
      energy%h = zero
      energy%h_old = zero
      energy%rhs_old = zero
      energy%rhs_old_valid = .false.
      energy%qrad = zero
      energy%cp = zero
      energy%lambda = zero
      energy%rho_thermo = zero
      energy%operator_consistent_rho_h = zero
      energy%operator_consistent_rho_h_available = 0
      energy%species_enthalpy_diffusion = zero
      energy%chem_heat_release_rate = zero
      energy%chem_dT_dt = zero
      energy%chem_dY_max_dt = zero
      energy%chem_rel_rho_change_dt = zero
      energy%chem_source_alpha = one
      energy%chem_active = zero
      energy%last_conductive_boundary_flux_out = zero
      energy%last_conductive_boundary_flux_available = 0
      energy%cumulative_boundary_rho_h_advective_flux_out = zero
      energy%cumulative_boundary_rho_h_conductive_flux_out = zero
      energy%cumulative_qrad_integral = zero
      energy%cumulative_rho_species_hdiff_integral = zero
      energy%cumulative_energy_budget_available = 0
      energy%last_energy_update_delta_rate_integral = zero
      energy%last_energy_update_rhs_integral = zero
      energy%last_energy_update_balance_defect = zero
      energy%relative_last_energy_update_balance_defect = zero
      energy%last_operator_consistent_rho_h_integral = zero
      energy%cumulative_energy_update_delta_integral = zero
      energy%cumulative_energy_update_rhs_integral = zero
      energy%initialized = .true.
   end subroutine allocate_energy