mod_energy.f90 Source File

!> Energy field storage, diagnostics, and sensible-enthalpy transport.
!!
!! This module owns the enthalpy/temperature/radiation-source fields used by
!! the current constant-density energy path. The transported thermodynamic
!! state is `h`; temperature is recovered from `h`, composition `Y`, and the
!! background thermodynamic pressure `params%background_press`.
!!
!! In Cantera mode the stored sensible enthalpy is
!! `h_abs(T,Y,p0) - h_abs(T_ref,Y,p0)`, where `p0` is the background pressure
!! and `T_ref` is `params%energy_reference_T`. After species transport changes
!! `Y`, Option A is enforced: preserve transported `h` and recover
!! `T(h,Y,p0)`. Do not rebuild `h` from the old temperature and new
!! composition.
!!
!! The energy equation uses constant flow/projection density `params%rho`.
!! `energy%rho_thermo` is diagnostic/future-use only and is not used by the
!! projection or momentum equations. Heat conduction uses `grad(T)`, not
!! `grad(h)`. `qrad` is a volumetric source with `qrad > 0` adding energy to
!! the gas.
!!
!! Missing physics in this path: reactions, reaction heat release,
!! variable-density low-Mach coupling, species-diffusion enthalpy correction
!! `-div(sum_k h_k J_k)`, and external radiation physics.
module mod_energy
   use mpi_f08
   use iso_c_binding, only : c_char, c_double, c_int
   use mod_kinds, only : rk, zero, tiny_safe, fatal_error, output_unit
   use mod_mesh_types, only : mesh_t
   use mod_input, only : case_params_t
   use mod_mpi_flow, only : flow_mpi_t, flow_exchange_cell_scalar
   use mod_bc, only : bc_set_t, bc_periodic, face_effective_neighbor, &
                      patch_type_for_face, boundary_temperature, boundary_species
   use mod_fields, only : flow_fields_t
   use mod_profiler, only : profiler_start, profiler_stop
   implicit none

   private

   !> Cell-centered energy variables.
   type, public :: energy_fields_t
      real(rk), allocatable :: T(:)          !< Temperature [K].
      real(rk), allocatable :: T_old(:)      !< Previous temperature state [K].
      real(rk), allocatable :: h(:)          !< Transported mixture sensible enthalpy [J/kg].
      real(rk), allocatable :: h_old(:)      !< Previous enthalpy state [J/kg].
      real(rk), allocatable :: qrad(:)       !< Volumetric source term [W/m^3]; positive values heat the gas.
      real(rk), allocatable :: cp(:)         !< Mixture heat capacity at constant pressure [J/kg/K].
      real(rk), allocatable :: lambda(:)     !< Thermal conductivity [W/m/K].
      real(rk), allocatable :: rho_thermo(:) !< Cantera thermodynamic density [kg/m^3], diagnostic/future-use.
      logical :: initialized = .false.  !< True after allocation/initialization.
   end type energy_fields_t

   public :: allocate_energy, initialize_energy, finalize_energy
   public :: update_enthalpy_from_temperature_constant_cp
   public :: recover_temperature_constant_cp
   public :: zero_radiation_source
   public :: advance_energy_transport
   public :: write_energy_diagnostics_header, write_energy_diagnostics_row

   !> C-Binding interface for Cantera thermodynamics.
   interface
      subroutine cantera_update_thermo_c(ncells, T, P, nspecies, Y_in, h_out, cp_out, lambda_out, &
                                         rho_thermo_out, T_ref, species_names_flat, name_len) bind(c, name="cantera_update_thermo_c")
         import :: c_char, c_double, c_int
         integer(c_int), value :: ncells
         real(c_double), intent(in) :: T(ncells), P(ncells)
         integer(c_int), value :: nspecies
         real(c_double), intent(in) :: Y_in(*)
         real(c_double), intent(out) :: h_out(ncells), cp_out(ncells), lambda_out(ncells), rho_thermo_out(ncells)
         real(c_double), value :: T_ref
         character(kind=c_char), intent(in) :: species_names_flat(*)
         integer(c_int), value :: name_len
      end subroutine cantera_update_thermo_c

      subroutine cantera_recover_temperature_from_h_c(ncells, h_in, P, nspecies, Y_in, T_out, &
                                                      T_ref, species_names_flat, name_len) bind(c, name="cantera_recover_temperature_from_h_c")
         import :: c_char, c_double, c_int
         integer(c_int), value :: ncells
         real(c_double), intent(in) :: h_in(ncells), P(ncells)
         integer(c_int), value :: nspecies
         real(c_double), intent(in) :: Y_in(*)
         real(c_double), intent(out) :: T_out(ncells)
         real(c_double), value :: T_ref
         character(kind=c_char), intent(in) :: species_names_flat(*)
         integer(c_int), value :: name_len
      end subroutine cantera_recover_temperature_from_h_c


      subroutine cantera_recover_temperature_and_update_thermo_c(ncells, h_in, P, nspecies, Y_in, T_out, &
                                                                 cp_out, lambda_out, rho_thermo_out, &
                                                                 T_ref, species_names_flat, name_len) &
         bind(c, name="cantera_recover_temperature_and_update_thermo_c")
         import :: c_char, c_double, c_int
         integer(c_int), value :: ncells
         real(c_double), intent(in) :: h_in(ncells), P(ncells)
         integer(c_int), value :: nspecies
         real(c_double), intent(in) :: Y_in(*)
         real(c_double), intent(out) :: T_out(ncells), cp_out(ncells), lambda_out(ncells), rho_thermo_out(ncells)
         real(c_double), value :: T_ref
         character(kind=c_char), intent(in) :: species_names_flat(*)
         integer(c_int), value :: name_len
      end subroutine cantera_recover_temperature_and_update_thermo_c

   end interface

contains

   !> Allocate all energy arrays for the mesh.
   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%qrad(mesh%ncells))
      allocate(energy%cp(mesh%ncells))
      allocate(energy%lambda(mesh%ncells))
      allocate(energy%rho_thermo(mesh%ncells))

      energy%T = zero
      energy%T_old = zero
      energy%h = zero
      energy%h_old = zero
      energy%qrad = zero
      energy%cp = zero
      energy%lambda = zero
      energy%rho_thermo = zero
      energy%initialized = .true.
   end subroutine allocate_energy


   !> Initialize energy fields from case parameters.
   subroutine initialize_energy(mesh, params, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      call allocate_energy(mesh, energy)

      energy%T = params%initial_T
      energy%T_old = energy%T

      if (params%enable_cantera_thermo) then
         call update_thermo_from_temperature_cantera(mesh, params, energy, species_Y)
      else
         call update_enthalpy_from_temperature_constant_cp(params, energy)
         energy%cp = params%energy_cp
         energy%lambda = params%energy_lambda
         energy%rho_thermo = params%rho
      end if

      energy%h_old = energy%h
      call zero_radiation_source(energy)
   end subroutine initialize_energy


   !> Deallocate all energy arrays.
   subroutine finalize_energy(energy)
      type(energy_fields_t), intent(inout) :: energy

      if (allocated(energy%T)) deallocate(energy%T)
      if (allocated(energy%T_old)) deallocate(energy%T_old)
      if (allocated(energy%h)) deallocate(energy%h)
      if (allocated(energy%h_old)) deallocate(energy%h_old)
      if (allocated(energy%qrad)) deallocate(energy%qrad)
      if (allocated(energy%cp)) deallocate(energy%cp)
      if (allocated(energy%lambda)) deallocate(energy%lambda)
      if (allocated(energy%rho_thermo)) deallocate(energy%rho_thermo)

      energy%initialized = .false.
   end subroutine finalize_energy


   !> Update h from T using the constant-cp thermodynamic model.
   subroutine update_enthalpy_from_temperature_constant_cp(params, energy)
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy

      if (.not. allocated(energy%T) .or. .not. allocated(energy%h)) then
         call fatal_error('energy', 'energy arrays are not allocated')
      end if

      if (params%energy_cp <= tiny_safe) then
         call fatal_error('energy', 'energy_cp must be positive')
      end if

      energy%h = params%energy_reference_h + &
                 params%energy_cp * (energy%T - params%energy_reference_T)

      if (allocated(energy%cp)) energy%cp = params%energy_cp
      if (allocated(energy%lambda)) energy%lambda = params%energy_lambda
      if (allocated(energy%rho_thermo)) energy%rho_thermo = params%rho
   end subroutine update_enthalpy_from_temperature_constant_cp


   !> Recover T from h using the constant-cp thermodynamic model.
   subroutine recover_temperature_constant_cp(params, energy)
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy

      if (.not. allocated(energy%T) .or. .not. allocated(energy%h)) then
         call fatal_error('energy', 'energy arrays are not allocated')
      end if

      if (params%energy_cp <= tiny_safe) then
         call fatal_error('energy', 'energy_cp must be positive')
      end if

      energy%T = params%energy_reference_T + &
                 (energy%h - params%energy_reference_h) / params%energy_cp
   end subroutine recover_temperature_constant_cp


   !> Constant-cp helper for a single boundary/face temperature value.
   pure function enthalpy_from_temperature_value(params, temperature) result(h_value)
      type(case_params_t), intent(in) :: params
      real(rk), intent(in) :: temperature
      real(rk) :: h_value

      h_value = params%energy_reference_h + &
                params%energy_cp * (temperature - params%energy_reference_T)
   end function enthalpy_from_temperature_value


   !> Build a flattened C-compatible species-name buffer.
   subroutine build_c_species_names(params, c_names_flat, n_len)
      type(case_params_t), intent(in) :: params
      character(kind=c_char), allocatable, intent(out) :: c_names_flat(:)
      integer, intent(out) :: n_len

      integer :: k, j, nsp

      nsp = max(1, params%nspecies)
      n_len = len(params%species_name(1))
      allocate(c_names_flat(n_len * nsp))
      c_names_flat = ' '

      do k = 1, nsp
         do j = 1, n_len
            c_names_flat((k-1)*n_len + j) = char(iachar(params%species_name(k)(j:j)), kind=c_char)
         end do
      end do
   end subroutine build_c_species_names


   !> Build a cellwise thermodynamic composition array for Cantera calls.
   !!
   !! If transported species are available, they are used directly and clipped
   !! to non-negative values. If their sum is positive, the local composition is
   !! normalized before passing to Cantera. If no transported species are
   !! available, use the single-species/default mixture prepared during Cantera
   !! initialization, usually thermo_default_species.
   subroutine build_thermo_Y(params, ncells, Y_local, species_Y)
      type(case_params_t), intent(in) :: params
      integer, intent(in) :: ncells
      real(rk), allocatable, intent(out) :: Y_local(:,:)
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: c, k, nsp
      real(rk) :: sum_Y

      nsp = max(1, params%nspecies)
      allocate(Y_local(nsp, ncells))
      Y_local = zero

      if (present(species_Y) .and. params%enable_species .and. params%nspecies > 0) then
         if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < ncells) then
            call fatal_error('energy', 'species_Y has incompatible shape for Cantera thermo update')
         end if

         do c = 1, ncells
            sum_Y = zero
            do k = 1, params%nspecies
               Y_local(k, c) = max(zero, species_Y(k, c))
               sum_Y = sum_Y + Y_local(k, c)
            end do

            if (sum_Y > tiny_safe) then
               Y_local(1:params%nspecies, c) = Y_local(1:params%nspecies, c) / sum_Y
            else
               Y_local(1, c) = 1.0_rk
            end if
         end do
      else
         Y_local(1, :) = 1.0_rk
      end if
   end subroutine build_thermo_Y


   !> Build a boundary thermodynamic composition vector for fixed-T boundaries.
   subroutine build_boundary_thermo_Y(mesh, bc, params, face_id, cell_id, species_Y, Y_point)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      real(rk), intent(in) :: species_Y(:,:)
      real(rk), intent(out) :: Y_point(:)

      integer :: k
      real(rk) :: ext_Y, sum_Y
      logical :: is_dirichlet

      if (size(Y_point) < params%nspecies) then
         call fatal_error('energy', 'Y_point too small in build_boundary_thermo_Y')
      end if
      if (size(species_Y, 1) < params%nspecies .or. size(species_Y, 2) < cell_id) then
         call fatal_error('energy', 'species_Y has incompatible shape in boundary thermo composition')
      end if

      sum_Y = zero
      do k = 1, params%nspecies
         call boundary_species(mesh, bc, face_id, k, species_Y(k, cell_id), ext_Y, is_dirichlet)
         Y_point(k) = max(zero, ext_Y)
         sum_Y = sum_Y + Y_point(k)
      end do

      if (sum_Y > tiny_safe) then
         Y_point(1:params%nspecies) = Y_point(1:params%nspecies) / sum_Y
      else
         Y_point = zero
         Y_point(1) = 1.0_rk
      end if
   end subroutine build_boundary_thermo_Y



   !> Update h, cp, lambda, and diagnostic thermo density from current T.
   subroutine update_thermo_from_temperature_cantera(mesh, params, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: n_len
      real(rk), allocatable :: P_arr(:), Y_local(:,:)
      character(kind=c_char), allocatable :: c_names_flat(:)

      if (.not. params%enable_cantera_thermo) return

      if (params%nspecies <= 0) then
         call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species')
      end if

      allocate(P_arr(mesh%ncells))
      P_arr = params%background_press
      call build_thermo_Y(params, mesh%ncells, Y_local, species_Y)
      call build_c_species_names(params, c_names_flat, n_len)

      call cantera_update_thermo_c(mesh%ncells, energy%T, P_arr, params%nspecies, Y_local, &
                                   energy%h, energy%cp, energy%lambda, energy%rho_thermo, &
                                   params%energy_reference_T, c_names_flat, n_len)

      deallocate(P_arr)
      deallocate(Y_local)
      deallocate(c_names_flat)
   end subroutine update_thermo_from_temperature_cantera


   !> Recover T from h using Cantera HPY inversion.
   subroutine recover_temperature_cantera(mesh, params, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: n_len
      real(rk), allocatable :: P_arr(:), Y_local(:,:)
      character(kind=c_char), allocatable :: c_names_flat(:)

      if (.not. params%enable_cantera_thermo) return

      if (params%nspecies <= 0) then
         call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species')
      end if

      allocate(P_arr(mesh%ncells))
      P_arr = params%background_press
      call build_thermo_Y(params, mesh%ncells, Y_local, species_Y)
      call build_c_species_names(params, c_names_flat, n_len)

      call cantera_recover_temperature_from_h_c(mesh%ncells, &
                                   energy%h, &
                                   P_arr, &
                                   params%nspecies, &
                                   Y_local, &
                                   energy%T, &
                                   params%energy_reference_T, &
                                   c_names_flat, &
                                   n_len)

      deallocate(P_arr)
      deallocate(Y_local)
      deallocate(c_names_flat)
   end subroutine recover_temperature_cantera


   !> Recover T from h and refresh cp/lambda/rho_thermo in one Cantera sync.
   !!
   !! This is the preferred energy-step thermo sync routine:
   !! `(T, cp, lambda, rho_thermo) = sync(h, Y, p0)`. It preserves the
   !! transported sensible enthalpy field and updates only derived thermo
   !! state: T, cp, lambda, and diagnostic rho_thermo.
   subroutine recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: n_len
      real(rk), allocatable :: P_arr(:), Y_local(:,:)
      character(kind=c_char), allocatable :: c_names_flat(:)

      if (.not. params%enable_cantera_thermo) return

      if (params%nspecies <= 0) then
         call fatal_error('energy', 'enable_cantera_thermo requires at least one inert thermo species')
      end if

      allocate(P_arr(mesh%ncells))
      P_arr = params%background_press
      call build_thermo_Y(params, mesh%ncells, Y_local, species_Y)
      call build_c_species_names(params, c_names_flat, n_len)

      call cantera_recover_temperature_and_update_thermo_c(mesh%ncells, &
                                   energy%h, &
                                   P_arr, &
                                   params%nspecies, &
                                   Y_local, &
                                   energy%T, &
                                   energy%cp, &
                                   energy%lambda, &
                                   energy%rho_thermo, &
                                   params%energy_reference_T, &
                                   c_names_flat, &
                                   n_len)

      deallocate(P_arr)
      deallocate(Y_local)
      deallocate(c_names_flat)
   end subroutine recover_temperature_and_update_thermo_cantera


   !> Compute h(T,Y,p0) for one boundary state using Cantera.
   subroutine enthalpy_from_temperature_cantera_point(params, temperature, h_value, Y_point)
      type(case_params_t), intent(in) :: params
      real(rk), intent(in) :: temperature
      real(rk), intent(out) :: h_value
      real(rk), intent(in), optional :: Y_point(:)

      integer :: n_len, nsp
      real(rk) :: T_arr(1), P_arr(1)
      real(rk) :: h_arr(1), cp_arr(1), lambda_arr(1), rho_arr(1)
      real(rk), allocatable :: Y_local(:,:)
      character(kind=c_char), allocatable :: c_names_flat(:)

      if (params%nspecies <= 0) then
         call fatal_error('energy', 'Cantera boundary enthalpy requires at least one species')
      end if

      nsp = max(1, params%nspecies)
      allocate(Y_local(nsp, 1))
      Y_local = zero

      if (present(Y_point)) then
         if (size(Y_point) < params%nspecies) then
            call fatal_error('energy', 'Y_point has incompatible size for Cantera boundary enthalpy')
         end if
         Y_local(1:params%nspecies, 1) = max(zero, Y_point(1:params%nspecies))
         if (sum(Y_local(1:params%nspecies, 1)) > tiny_safe) then
            Y_local(1:params%nspecies, 1) = Y_local(1:params%nspecies, 1) / &
                                            sum(Y_local(1:params%nspecies, 1))
         else
            Y_local(1, 1) = 1.0_rk
         end if
      else
         Y_local(1, 1) = 1.0_rk
      end if

      T_arr(1) = temperature
      P_arr(1) = params%background_press
      call build_c_species_names(params, c_names_flat, n_len)

      call cantera_update_thermo_c(1, T_arr, P_arr, params%nspecies, Y_local, &
                                   h_arr, cp_arr, lambda_arr, rho_arr, &
                                   params%energy_reference_T, c_names_flat, n_len)

      h_value = h_arr(1)

      deallocate(Y_local)
      deallocate(c_names_flat)
   end subroutine enthalpy_from_temperature_cantera_point



   !> Convert a boundary temperature to enthalpy using the active thermo model.
   
   !> Boundary enthalpy from a boundary temperature.
   !!
   !! When a transported boundary composition is available, the fixed-temperature
   !! boundary enthalpy is evaluated as h(T_b,Y_b,p0). This keeps hot/cold
   !! species inlets thermodynamically consistent.
   function boundary_enthalpy_from_temperature(mesh, params, temperature, Y_point) result(h_value)
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params
      real(rk), intent(in) :: temperature
      real(rk), intent(in), optional :: Y_point(:)
      real(rk) :: h_value

      if (.not. params%enable_cantera_thermo) then
         h_value = enthalpy_from_temperature_value(params, temperature)
         return
      end if

      call enthalpy_from_temperature_cantera_point(params, temperature, h_value, Y_point)

      associate(dummy => mesh%ncells)
      end associate
   end function boundary_enthalpy_from_temperature



   !> Reset the volumetric energy source to zero.
   subroutine zero_radiation_source(energy)
      type(energy_fields_t), intent(inout) :: energy

      if (allocated(energy%qrad)) energy%qrad = zero
   end subroutine zero_radiation_source


   !> Advance transported sensible enthalpy with constant flow density.
   !!
   !! The explicit finite-volume update is:
   !!
   !!   V dh/dt = - sum_f F_f h_f
   !!            + (1/rho) sum_f lambda A_f (T_nb - T_c)/d_f
   !!            + (qrad/rho) V
   !!
   !! fields%face_flux is volumetric flux. It is re-oriented outward from the
   !! currently updated cell before applying upwind h advection. Cantera thermo
   !! may provide `cp` and `lambda`, but the projection/momentum density remains
   !! `params%rho`.
   subroutine advance_energy_transport(mesh, flow, bc, params, fields, energy, species_Y)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(inout) :: energy
      real(rk), intent(in), optional :: species_Y(:,:)

      integer :: c, lf, fid, neighbor
      real(rk) :: flux_out
      real(rk) :: face_area, dist
      real(rk) :: h_cell, h_other, h_face
      real(rk) :: T_cell, T_other
      real(rk) :: rhs_h, diff_term, lambda_face
      logical :: is_dirichlet, do_diffusion
      real(rk), allocatable :: Y_point(:)

      if (.not. params%enable_energy) return

      if (.not. energy%initialized) then
         call fatal_error('energy', 'advance requested before energy initialization')
      end if

      if (params%rho <= tiny_safe) call fatal_error('energy', 'rho must be positive')
      if (params%energy_cp <= tiny_safe) call fatal_error('energy', 'energy_cp must be positive')
      if (params%energy_lambda < zero) call fatal_error('energy', 'energy_lambda must be non-negative')

      ! Option A: preserve transported enthalpy across species updates.
      ! Species transport may have changed Y before this energy step.  Do not
      ! recompute h from the old temperature and the new composition.  Instead,
      ! keep h as the transported state, recover T(h,Y,p0), then refresh the
      ! thermo/transport properties used by the enthalpy equation.
      call profiler_start('Energy_Exchange_H')
      call flow_exchange_cell_scalar(flow, energy%h)
      call profiler_stop('Energy_Exchange_H')

      energy%h_old = energy%h

      if (params%enable_cantera_thermo) then
         ! A pre-flux thermo sync is required only when species transport may
         ! have changed the composition since the previous post-flux sync.
         ! Without species, current T/cp/lambda/rho_thermo are already valid
         ! from initialization or the previous energy step's post-sync.
         if (params%enable_species .and. present(species_Y) .and. params%nspecies > 0) then
            call profiler_start('Energy_Cantera_PreSync')
            call recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y)
            call profiler_stop('Energy_Cantera_PreSync')

            energy%h = energy%h_old
         end if
      end if

      ! Ensure off-rank temperature/property values are current before fluxes.
      call profiler_start('Energy_PreFlux_Exchange')
      call flow_exchange_cell_scalar(flow, energy%T)
      if (allocated(energy%lambda)) call flow_exchange_cell_scalar(flow, energy%lambda)
      call profiler_stop('Energy_PreFlux_Exchange')

      if (params%enable_cantera_thermo .and. params%enable_species .and. &
          present(species_Y) .and. params%nspecies > 0) then
         allocate(Y_point(params%nspecies))
      end if
      if (allocated(energy%T_old)) energy%T_old = energy%T

      call profiler_start('Energy_Flux_Update')
      do c = flow%first_cell, flow%last_cell
         rhs_h = zero
         h_cell = energy%h_old(c)
         T_cell = energy%T(c)

         do lf = 1, mesh%ncell_faces(c)
            fid = mesh%cell_faces(lf, c)

            ! fields%face_flux is oriented owner -> neighbor. Reorient outward
            ! from the current cell.
            if (mesh%faces(fid)%owner == c) then
               flux_out = fields%face_flux(fid)
            else
               flux_out = -fields%face_flux(fid)
            end if

            neighbor = face_effective_neighbor(mesh, bc, fid, c)
            face_area = mesh%faces(fid)%area

            if (neighbor > 0) then
               h_other = energy%h_old(neighbor)
               T_other = energy%T(neighbor)
               do_diffusion = .true.
            else
               call boundary_temperature(mesh, bc, fid, T_cell, T_other, is_dirichlet)
               if (is_dirichlet) then
                  if (allocated(Y_point)) then
                     call build_boundary_thermo_Y(mesh, bc, params, fid, c, species_Y, Y_point)
                     h_other = boundary_enthalpy_from_temperature(mesh, params, T_other, Y_point)
                  else
                     h_other = boundary_enthalpy_from_temperature(mesh, params, T_other)
                  end if
                  do_diffusion = .true.
               else
                  h_other = h_cell
                  T_other = T_cell
                  do_diffusion = .false.
               end if
            end if

            ! Upwind advection of h using outward volumetric flux.
            if (flux_out >= zero) then
               h_face = h_cell
            else
               h_face = h_other
            end if
            rhs_h = rhs_h - flux_out * h_face

            ! Fourier heat conduction uses grad(T), not grad(h).
            if (do_diffusion .and. (params%enable_cantera_thermo .or. params%energy_lambda > zero)) then
               dist = energy_face_normal_distance(mesh, bc, fid, c, neighbor)

               if (allocated(energy%lambda)) then
                  if (neighbor > 0) then
                     lambda_face = 0.5_rk * (energy%lambda(c) + energy%lambda(neighbor))
                  else
                     lambda_face = energy%lambda(c)
                  end if
               else
                  lambda_face = params%energy_lambda
               end if

               if (lambda_face > zero) then
                  diff_term = (lambda_face / params%rho) * &
                              (T_other - T_cell) / dist * face_area
                  rhs_h = rhs_h + diff_term
               end if
            end if
         end do

         energy%h(c) = energy%h_old(c) + params%dt * &
                       (rhs_h / mesh%cells(c)%volume + energy%qrad(c) / params%rho)
      end do
      call profiler_stop('Energy_Flux_Update')

      if (params%enable_cantera_thermo) then
         ! Protect transported h from property-refresh roundoff: recover T from
         ! the newly transported h, refresh cp/lambda/rho_thermo, then restore h.
         energy%h_old = energy%h

         call profiler_start('Energy_Cantera_PostSync')
         call recover_temperature_and_update_thermo_cantera(mesh, params, energy, species_Y)
         call profiler_stop('Energy_Cantera_PostSync')

         energy%h = energy%h_old
      else
         call profiler_start('Energy_ConstantCp_RecoverT')
         call recover_temperature_constant_cp(params, energy)
         call profiler_stop('Energy_ConstantCp_RecoverT')
      end if

      ! Synchronize updated owned cells for output and the next step.
      call profiler_start('Energy_Final_Exchange')
      call flow_exchange_cell_scalar(flow, energy%h)
      call flow_exchange_cell_scalar(flow, energy%T)
      if (allocated(energy%lambda)) call flow_exchange_cell_scalar(flow, energy%lambda)
      call profiler_stop('Energy_Final_Exchange')
   end subroutine advance_energy_transport


   !> Outward unit normal from cell_id for face_id.
   function energy_outward_normal(mesh, face_id, cell_id) result(nvec)
      type(mesh_t), intent(in) :: mesh
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      real(rk) :: nvec(3)

      if (mesh%faces(face_id)%owner == cell_id) then
         nvec = mesh%faces(face_id)%normal
      else
         nvec = -mesh%faces(face_id)%normal
      end if
   end function energy_outward_normal


   !> Normal distance used for cell-cell and cell-boundary temperature gradients.
   function energy_face_normal_distance(mesh, bc, face_id, cell_id, nb) result(dist)
      type(mesh_t), intent(in) :: mesh
      type(bc_set_t), intent(in) :: bc
      integer, intent(in) :: face_id
      integer, intent(in) :: cell_id
      integer, intent(in) :: nb
      real(rk) :: dist

      integer :: pair_face
      integer :: btype
      real(rk) :: nvec(3)

      nvec = energy_outward_normal(mesh, face_id, cell_id)

      if (nb > 0) then
         if (mesh%faces(face_id)%neighbor == 0) then
            btype = patch_type_for_face(mesh, bc, face_id)

            if (btype == bc_periodic) then
               pair_face = mesh%faces(face_id)%periodic_face

               if (pair_face <= 0) then
                  call fatal_error('energy', 'periodic face has no paired face')
               end if

               dist = abs(dot_product(mesh%faces(face_id)%center - &
                                      mesh%cells(cell_id)%center, nvec)) + &
                      abs(dot_product(mesh%cells(nb)%center - &
                                      mesh%faces(pair_face)%center, nvec))

               dist = max(dist, tiny_safe)
               return
            end if
         end if

         dist = abs(dot_product(mesh%cells(nb)%center - &
                                mesh%cells(cell_id)%center, nvec))
      else
         dist = abs(dot_product(mesh%faces(face_id)%center - &
                                mesh%cells(cell_id)%center, nvec))
      end if

      dist = max(dist, tiny_safe)
   end function energy_face_normal_distance


   !> Writes the CSV header for energy diagnostics.
   subroutine write_energy_diagnostics_header(params, flow)
      type(case_params_t), intent(in) :: params
      type(flow_mpi_t), intent(in) :: flow

      integer :: unit_id
      character(len=256 + 32) :: filename

      if (flow%rank /= 0 .or. .not. params%write_diagnostics) return
      if (.not. params%enable_energy) return

      filename = trim(params%output_dir)//'/energy_diagnostics.csv'
      open(newunit=unit_id, file=trim(filename), status='replace', action='write')

      write(unit_id,'(a)') 'step,time,min_T,max_T,mean_T,min_h,max_h,mean_h,min_qrad,max_qrad,integral_qrad,max_delta_T,rel_h_residual'

      close(unit_id)
   end subroutine write_energy_diagnostics_header


   !> Appends one row of global energy diagnostics.
   !!
   !! The temperature and enthalpy means are volume weighted over owned cells.
   !! qrad integral is the domain integral of qrad dV [W].
   subroutine write_energy_diagnostics_row(mesh, flow, params, energy, step, time)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(energy_fields_t), intent(in) :: energy
      integer, intent(in) :: step
      real(rk), intent(in) :: time

      integer :: c, unit_id, ierr
      character(len=256 + 32) :: filename
      real(rk) :: vol, old_T
      real(rk) :: local_min_T, local_max_T, local_sum_T
      real(rk) :: local_min_h, local_max_h, local_sum_h
      real(rk) :: local_min_qrad, local_max_qrad, local_integral_qrad
      real(rk) :: local_volume
      real(rk) :: local_max_delta_T, local_delta_h2, local_h_old2
      real(rk) :: global_min_T, global_max_T, global_sum_T
      real(rk) :: global_min_h, global_max_h, global_sum_h
      real(rk) :: global_min_qrad, global_max_qrad, global_integral_qrad
      real(rk) :: global_volume
      real(rk) :: global_max_delta_T, global_delta_h2, global_h_old2
      real(rk) :: mean_T, mean_h, rel_h_residual
      logical, save :: printed_thermo_mode = .false.
      if (.not. params%write_diagnostics) return
      if (.not. params%enable_energy) return

      if (.not. allocated(energy%T) .or. .not. allocated(energy%h) .or. &
          .not. allocated(energy%h_old) .or. .not. allocated(energy%qrad)) then
         call fatal_error('energy', 'energy diagnostics requested but energy arrays are not allocated')
      end if

      local_min_T = huge(0.0_rk)
      local_max_T = -huge(0.0_rk)
      local_sum_T = zero
      local_min_h = huge(0.0_rk)
      local_max_h = -huge(0.0_rk)
      local_sum_h = zero
      local_min_qrad = huge(0.0_rk)
      local_max_qrad = -huge(0.0_rk)
      local_integral_qrad = zero
      local_volume = zero
      local_max_delta_T = zero
      local_delta_h2 = zero
      local_h_old2 = zero

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

         vol = mesh%cells(c)%volume
         local_volume = local_volume + vol

         local_min_T = min(local_min_T, energy%T(c))
         local_max_T = max(local_max_T, energy%T(c))
         local_sum_T = local_sum_T + energy%T(c) * vol

         local_min_h = min(local_min_h, energy%h(c))
         local_max_h = max(local_max_h, energy%h(c))
         local_sum_h = local_sum_h + energy%h(c) * vol

         local_min_qrad = min(local_min_qrad, energy%qrad(c))
         local_max_qrad = max(local_max_qrad, energy%qrad(c))
         local_integral_qrad = local_integral_qrad + energy%qrad(c) * vol

         if (allocated(energy%T_old)) then
            old_T = energy%T_old(c)
         else
            old_T = params%energy_reference_T + &
                    (energy%h_old(c) - params%energy_reference_h) / params%energy_cp
         end if
         local_max_delta_T = max(local_max_delta_T, abs(energy%T(c) - old_T))
         local_delta_h2 = local_delta_h2 + (energy%h(c) - energy%h_old(c))**2 * vol
         local_h_old2 = local_h_old2 + energy%h_old(c)**2 * vol
      end do

      call MPI_Allreduce(local_min_T, global_min_T, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_T')
      call MPI_Allreduce(local_max_T, global_max_T, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_T')
      call MPI_Allreduce(local_sum_T, global_sum_T, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing sum_T')

      call MPI_Allreduce(local_min_h, global_min_h, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_h')
      call MPI_Allreduce(local_max_h, global_max_h, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_h')
      call MPI_Allreduce(local_sum_h, global_sum_h, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing sum_h')

      call MPI_Allreduce(local_min_qrad, global_min_qrad, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing min_qrad')
      call MPI_Allreduce(local_max_qrad, global_max_qrad, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_qrad')
      call MPI_Allreduce(local_integral_qrad, global_integral_qrad, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing integral_qrad')

      call MPI_Allreduce(local_volume, global_volume, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing energy volume')
      call MPI_Allreduce(local_max_delta_T, global_max_delta_T, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing max_delta_T')
      call MPI_Allreduce(local_delta_h2, global_delta_h2, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing delta_h2')
      call MPI_Allreduce(local_h_old2, global_h_old2, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('energy', 'MPI failure reducing h_old2')

      if (global_volume > tiny_safe) then
         mean_T = global_sum_T / global_volume
         mean_h = global_sum_h / global_volume
      else
         mean_T = zero
         mean_h = zero
      end if

      rel_h_residual = sqrt(global_delta_h2) / max(sqrt(global_h_old2), tiny_safe)

      if (flow%rank /= 0) return

      filename = trim(params%output_dir)//'/energy_diagnostics.csv'
      open(newunit=unit_id, file=trim(filename), status='old', position='append', action='write')

      write(unit_id,'(i0,12(",",es16.8))') step, time, global_min_T, global_max_T, mean_T, &
         global_min_h, global_max_h, mean_h, global_min_qrad, global_max_qrad, &
         global_integral_qrad, global_max_delta_T, rel_h_residual

      close(unit_id)
      if (params%enable_cantera_thermo .and. step > 0 .and. .not. printed_thermo_mode) then
         if (params%enable_species) then
            write(output_unit,'(a,i0)') 'Cantera thermo mode: transported species composition, nspecies = ', params%nspecies
         else
            write(output_unit,'(a,a)') 'Cantera thermo mode: default single species = ', trim(params%thermo_default_species)
         end if
         printed_thermo_mode = .true.
      end if

      write(output_unit,'(a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5)') &
         'energy step=', step, 'Tmin=', global_min_T, 'Tmax=', global_max_T, &
         'Tmean=', mean_T, 'dTmax=', global_max_delta_T, 'rel_h=', rel_h_residual
   end subroutine write_energy_diagnostics_row

end module mod_energy