!> Main entry point for the LowMachReact-Hex finite-volume solver.
!!
!! This program orchestrates the entire simulation lifecycle:
!! 1. **Initialization**: Starts MPI, parses the `case.nml` namelist, and reads the mesh.
!! 2. **Domain Decomposition**: Sets up the replicated mesh MPI ranks for flow and radiation.
!! 3. **Field Setup**: Allocates and initializes velocity, pressure, species, and energy fields.
!! 4. **Time Integration Loop**: Executes projection, non-reacting species transport,
!!    and sensible-enthalpy transport. Constant-density mode is the stable
!!    baseline; guarded variable-density mode syncs active density from
!!    Cantera thermo density and projects to \( \nabla \cdot u=S \).
!! 5. **Diagnostics & Output**: Computes global observables and writes VTU/CSV files.
!! 6. **Finalization**: Safely releases all allocated memory and shuts down MPI.
program lowmach_react_hex
   use mpi_f08
   use mod_precision, only : rk, zero, one, tiny_safe, output_unit, fatal_error
   use mod_input, only : case_params_t, read_case_params, max_species
   use mod_mesh, only : mesh_t, mesh_finalize
   use mod_mesh_io, only : read_native_mesh
   use mod_mpi_flow, only : flow_mpi_t, mpi_flow_startup, mpi_flow_shutdown, &
                            flow_mpi_initialize, flow_mpi_finalize, &
                            flow_allgather_owned_scalar, flow_allgather_owned_matrix
   use mod_mpi_radiation, only : radiation_mpi_t, radiation_mpi_initialize, &
                                 radiation_mpi_finalize, radiation_configure_spectral
   use mod_bc, only : bc_set_t, build_bc_set, finalize_bc_set
   use mod_species_composition, only : parse_named_composition, composition_string_is_set
   use mod_flow_fields, only : flow_fields_t, initialize_fields, finalize_fields
   use mod_flow_projection, only : solver_stats_t, advance_projection_step, compute_flow_diagnostics, compute_and_update_cfl, &
                                   update_lowmach_divergence_source_from_density, correct_projection_to_lowmach_source
   use mod_output, only : &
                           prepare_output, &
                           write_diagnostics_header, &
                           write_diagnostics_row, &
                           write_vtu_unstructured, &
                           write_mesh_summary, &
                           write_pvd_collection, &
                           write_variable_density_diagnostics, &
                           write_boundary_mass_balance_diagnostics, &
                           write_species_energy_conservation_diagnostics, &
                           load_existing_pvd
   use mod_species, only : species_fields_t, initialize_species, finalize_species, &
                           advance_species_transport
   use mod_energy, only : &
                           energy_fields_t, &
                           initialize_energy, &
                           finalize_energy, &
                           refresh_energy_from_temperature, &
                           advance_energy_transport, &
                           advance_cantera_chemistry, &
                           write_cantera_cache_stats, write_chemistry_screening_stats, &
                           write_energy_diagnostics_header, &
                           write_energy_diagnostics_row, &
                           write_chemistry_diagnostics_header
   use mod_cantera, only : transport_properties_t, initialize_transport, &
                                        finalize_transport, update_transport_properties, sync_active_density_from_thermo
   use mod_restart, only : read_restart_file, write_restart_file
   use mod_radiation, only : radiation_context_t, radiation_initialize, radiation_update_source, radiation_finalize
   use mod_profiling, only : profiler_start, profiler_stop, profiler_report, profiler_configure, &
                          profiler_register_standard_timers
   implicit none

   type :: step_backup_t
      type(flow_fields_t) :: fields
      type(species_fields_t) :: species
      type(energy_fields_t) :: energy
      type(transport_properties_t) :: transport
      type(solver_stats_t) :: stats
      real(rk) :: time = zero
      real(rk) :: dt = zero
      real(rk) :: chemistry_accum_dt = zero
      logical :: valid = .false.
   end type step_backup_t

   type(case_params_t) :: params    !< Parsed simulation parameters.
   type(mesh_t) :: mesh             !< Global computational mesh.
   type(flow_mpi_t) :: flow_mpi     !< MPI context for flow solver.
   type(radiation_mpi_t) :: rad_mpi !< MPI context for radiation solver.
   type(radiation_context_t) :: rad_context !< Radiation model/data-exchange context.
   type(bc_set_t) :: bc             !< Boundary condition definitions.
   type(flow_fields_t) :: fields    !< Hydrodynamic fields (U, P).
   type(solver_stats_t) :: stats    !< Runtime diagnostics and solver performance metrics.
   type(species_fields_t) :: species !< Species mass fractions.
   type(energy_fields_t) :: energy   !< Enthalpy/temperature/radiation-source fields.
   type(transport_properties_t) :: transport !< Physical properties (rho, mu, Dk).

   character(len=256) :: case_file
   integer :: step, start_step, ierr, world_rank
   logical :: do_output_step, do_terminal_step, do_status_step, do_cfl_step
   logical :: do_terminal_full, do_terminal_health
   logical :: is_chem_step, is_rad_step
   real(rk) :: time, t0, chemistry_accum_dt

   real(rk) :: console_max_divu, console_sproj_max
   real(rk) :: console_lm_res_projection, console_lm_res_current
   real(rk) :: console_mass_balance_defect
   real(rk) :: terminal_projection_metric, terminal_ke, terminal_tmin, terminal_tmax
   real(rk) :: chem_max_dT, chem_max_dY, chem_max_rel_drho, chem_min_alpha
   real(rk) :: chem_raw_max_dT, chem_raw_max_dY, chem_raw_max_rel_drho
   type(step_backup_t) :: step_backup
   integer :: step_retry_count
   logical :: step_accepted, reject_step, block_dynamic_dt
   character(len=256) :: reject_reason
   real(rk) :: step_ke_before, step_ke_after, step_cfl_after, step_umax_after
   real(rk) :: step_trial_dt, step_cut_dt
   ! 1. Initialize MPI and read runtime configuration.
   call mpi_flow_startup()
   call MPI_Comm_rank(MPI_COMM_WORLD, world_rank, ierr)
   if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI_Comm_rank failed')

   t0 = real(mpi_wtime(), rk)
   step = 0
   start_step = 0
   time = zero
   chemistry_accum_dt = zero

   call get_case_filename(case_file)
   call read_case_params(trim(case_file), params)
   call profiler_configure(params%enable_profiling, params%nested_profiling)
   call profiler_register_standard_timers()

   if (world_rank == 0) then
      write(output_unit,'(a,l1,a,l1,a,l1)') 'profiling config: enable=', &
         params%enable_profiling, ' nested=', params%nested_profiling, &
         ' cantera_cache_stats=', params%write_cantera_cache_stats
   end if

   call profiler_start('Total_Simulation')

   ! 2. Setup simulation environment and read input data.
   call read_native_mesh(trim(params%mesh_dir), mesh)
   
   ! 3. Initialize parallel contexts.
   call initialize_transport(mesh, params, transport)
   call flow_mpi_initialize(mesh, flow_mpi, MPI_COMM_WORLD, max(4, params%nspecies))
   call radiation_mpi_initialize(rad_mpi, MPI_COMM_WORLD)
   if (params%enable_radiation) then
      call radiation_configure_spectral(rad_mpi, params%radiation_n_wavenumbers)
   end if
   
   ! 4. Prepare physical fields and BCs.
   call build_bc_set(mesh, params, bc)
   call initialize_fields(mesh, fields)
   if (params%enable_species) then
      call initialize_species(mesh, params, species)
   end if
   if (params%enable_energy) then
      if (params%enable_species) then
         call initialize_energy(mesh, params, energy, species%Y)
      else
         call initialize_energy(mesh, params, energy)
      end if

      if (params%enable_ignition_kernel .and. .not. params%restart_from_file) then
         call apply_ignition_kernel(mesh, flow_mpi, params, species, energy)
         if (params%enable_species) then
            call refresh_energy_from_temperature(mesh, params, energy, species%Y)
         else
            call refresh_energy_from_temperature(mesh, params, energy)
         end if
      end if

      ! No-op for the default constant-density solver. In guarded
      ! density_eos="cantera" variable-density mode, energy%rho_thermo becomes
      ! the active flow density after thermo sync.
      call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo)

      ! Initial variable-density state: rho_old must equal active rho at t=0.
      ! The sync routine stores the pre-sync placeholder density in rho_old;
      ! that is useful during timestepping but not for the initial state.
      if (params%enable_variable_density) then
         if (allocated(transport%rho_old)) transport%rho_old = transport%rho
      end if
   end if

   ! Optional rank-independent restart.  The file stores global arrays in mesh
   ! order, so this works even when the new run uses a different MPI rank count.
   call read_restart_file(mesh, flow_mpi, params, fields, species, energy, transport, start_step, time, chemistry_accum_dt)
    
   ! If ignition kernel is requested on a restarted run, overwrite the fields
   ! immediately after loading and re-refresh enthalpy.
   if (params%enable_ignition_kernel .and. params%restart_from_file) then
      call apply_ignition_kernel(mesh, flow_mpi, params, species, energy)
      if (params%enable_species) then
         call refresh_energy_from_temperature(mesh, params, energy, species%Y)
      else
         call refresh_energy_from_temperature(mesh, params, energy)
      end if
   end if
    
   ! Align active flow density with thermodynamic density immediately after restart
   ! to prevent catastrophic startup history-source transients.
   if (params%enable_variable_density) then
      call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo)
      if (allocated(transport%rho_old)) transport%rho_old = transport%rho
   end if

   params%dt_old = params%dt

   if (flow_mpi%rank == 0 .and. params%enable_species) then
      print *, "species: ", species%nspecies
   end if

   ! 5. Setup output files.
   call prepare_output(params, flow_mpi)
   call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/diagnostics')
   call execute_command_line('mkdir -p ' // trim(params%output_dir) // '/VTK')
   call write_mesh_summary(params, flow_mpi, mesh)
   call write_diagnostics_header(params, flow_mpi)
   call write_energy_diagnostics_header(params, flow_mpi)
   call write_chemistry_diagnostics_header(params, flow_mpi)
   call load_existing_pvd(params, flow_mpi, start_step)
   call radiation_initialize(mesh, rad_mpi, params, rad_context)
   call write_run_config_summary(mesh, flow_mpi, rad_mpi, params, trim(case_file))

   ! 6. Initial diagnostics and visualization snapshot.
   step = start_step
   call compute_flow_diagnostics(mesh, flow_mpi, bc, params, fields, stats)
   call write_diagnostics_row(params, flow_mpi, step, time, stats)
   call write_boundary_mass_balance_diagnostics(mesh, flow_mpi, params, fields, transport, step, time)
   if (params%enable_variable_density .and. step > 0 .and. time > 0.0_rk) then
      call write_variable_density_diagnostics(mesh, flow_mpi, params, fields, energy, transport, step, time)
   end if
   call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, step, time, write_file=.true.)
   call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, step, time)

   if (flow_mpi%rank == 0) then
      write(output_unit,'(a)') 'LowMachReact-Hex hexahedral low-Mach FV solver'
      write(output_unit,'(a,a)') 'case: ', trim(case_file)
      write(output_unit,'(a,i0)') 'cells: ', mesh%ncells
      write(output_unit,'(a,i0)') 'flow/chemistry MPI ranks: ', flow_mpi%nprocs
      write(output_unit,'(a,i0)') 'radiation MPI ranks: ', rad_mpi%nprocs
      if (params%enable_radiation) then
         write(output_unit,'(a,i0,a,i0,a,i0)') 'Radiation spectral tasks: n_wavenumbers=', &
            params%radiation_n_wavenumbers, ' local(first)=', rad_mpi%first_wavenumber, &
            ' update_interval=', params%radiation_update_interval
         write(output_unit,'(a,a)') 'Radiation source model: ', trim(params%radiation_source_model)
         write(output_unit,'(a)') 'Radiation coupling: spectral communicator computes qrad on full replicated mesh and reuses qrad between updates'
      end if
      if (params%enable_variable_density) then
         write(output_unit,'(a,a)') 'Flow density mode: EXPERIMENTAL variable-density low-Mach, density_eos = ', &
            trim(params%density_eos)
         write(output_unit,'(a)') 'Flow density source: transport%rho <- energy%rho_thermo after thermo sync'
         write(output_unit,'(a)') 'Low-Mach constraint: div(u) = S with conservative density-history plus advective-density source'
         if (params%variable_density_debug) then
            write(output_unit,'(a)') 'Variable-density debug diagnostics: enabled'
         else
            write(output_unit,'(a)') 'Variable-density debug diagnostics: disabled'
         end if
      else
         write(output_unit,'(a,es12.5,a)') 'Flow density mode: constant rho = ', params%rho, ' kg/m^3'
      end if
      if (params%enable_variable_nu) then
         write(output_unit,'(a)') 'Flow viscosity mode: variable Cantera mu/nu enabled'
      else
         write(output_unit,'(a,es12.5,a)') 'Flow viscosity mode: constant nu = ', params%nu, ' m^2/s'
      end if
      if (params%enable_cantera_fluid .or. params%enable_cantera_species .or. params%enable_cantera_thermo) then
         if (len_trim(params%cantera_phase_name) > 0) then
            write(output_unit,'(a,a)') 'Cantera phase selector: ', trim(params%cantera_phase_name)
         else
            write(output_unit,'(a)') 'Cantera phase selector: <default first phase>'
         end if
      end if
      if (params%enable_energy .and. params%enable_cantera_thermo) then
         if (params%enable_variable_density .and. trim(params%density_eos) == 'cantera') then
            write(output_unit,'(a)') 'Cantera rho_thermo: ACTIVE EOS density for experimental low-Mach coupling'
         else
            write(output_unit,'(a)') 'Cantera rho_thermo: diagnostic only, not used by projection'
         end if
      end if
      if (params%enable_cantera_fluid .or. params%enable_cantera_species) then
         write(output_unit,'(a,i0,a)') 'Cantera transport update interval: ', &
            params%transport_update_interval, ' step(s) [mu/D_k only]'
         if (params%enable_energy) then
            write(output_unit,'(a)') 'Cantera transport temperature source: energy%T'
         else
            write(output_unit,'(a,es12.5,a)') 'Cantera transport temperature source: background_temp = ', &
               params%background_temp, ' K'
         end if
      end if
      if (params%enable_energy .and. params%enable_species_enthalpy_diffusion) then
         if (params%enable_variable_density) then
            write(output_unit,'(a)') 'Species enthalpy diffusion: enabled with density-weighted variable-density branch'
         else
            write(output_unit,'(a)') 'Species enthalpy diffusion: enabled [-div(sum_k h_k J_k)/rho]'
         end if
      end if
      if (params%enable_energy .and. params%enable_cantera_thermo) then
         write(output_unit,'(a)') 'Cantera thermo update interval: every energy step'
         write(output_unit,'(a,es12.5,a)') 'Cantera thermodynamic pressure p0: ', &
            params%background_press, ' Pa'
      end if
      if (params%enable_reactions) then
         write(output_unit,'(a,i0,a,es10.3,a,es10.3,a,i0)') 'Cantera chemistry: enabled, update_interval=', &
            params%chemistry_update_interval, ' rtol=', params%chemistry_rtol, &
            ' atol=', params%chemistry_atol, ' max_steps=', params%chemistry_max_steps
      end if
   end if

   ! 7. Main Time-Stepping Loop
   do step = start_step + 1, params%nsteps
      step_retry_count = 0
      step_accepted = .false.

      do while (.not. step_accepted)
         if (params%enable_step_rejection) then
            call save_step_backup(fields, species, energy, transport, stats, time, params%dt, chemistry_accum_dt, step_backup)
         end if
         call compute_step_health(mesh, flow_mpi, params, fields, transport, step_ke_before, step_cfl_after, step_umax_after)
         chem_max_dT = zero
         chem_max_dY = zero
         chem_max_rel_drho = zero
         chem_min_alpha = 1.0_rk
         chem_raw_max_dT = zero
         chem_raw_max_dY = zero
         chem_raw_max_rel_drho = zero

      do_output_step = (mod(step, params%output_interval) == 0) .or. step == params%nsteps
      do_terminal_step = (mod(step, params%terminal_interval) == 0) .or. step == params%nsteps
      do_terminal_full = trim(params%terminal_detail) == 'full'
      do_terminal_health = trim(params%terminal_detail) == 'health'
      do_status_step = do_output_step .or. (do_terminal_step .and. do_terminal_full)
      do_cfl_step = params%use_dynamic_dt .or. do_output_step .or. &
                    (do_terminal_step .and. (do_terminal_health .or. do_terminal_full))

      ! A. Adaptive time-stepping / CFL diagnostics.
      ! Dynamic dt requires CFL every step. Health/full terminal modes also
      ! report CFL, while brief mode avoids this diagnostic reduction unless
      ! output or dynamic timestep control needs it.
      if (do_cfl_step) then
         call profiler_start('CFL_Update')
         if (step_retry_count > 0) then
            block_dynamic_dt = params%use_dynamic_dt
            params%use_dynamic_dt = .false.
            call compute_and_update_cfl(mesh, flow_mpi, params, fields, stats)
            params%use_dynamic_dt = block_dynamic_dt
         else
            call compute_and_update_cfl(mesh, flow_mpi, params, fields, stats)
         end if
         call profiler_stop('CFL_Update')
      end if

      ! B. Dynamic transport-property update via Cantera.
      ! transport_update_interval gates Cantera mu/D_k refresh; enable_variable_nu controls whether mu affects flow.
      ! In variable-density mode, active flow density is refreshed after the energy thermo sync.
      if (mod(step-1, params%transport_update_interval) == 0 .or. step == 1) then
         call profiler_start('Transport_Update')
         if (params%enable_species) then
            if (params%enable_energy) then
               call update_transport_properties(mesh, flow_mpi, params, species%Y, transport, T_state=energy%T)
            else
               call update_transport_properties(mesh, flow_mpi, params, species%Y, transport)
            end if
         else
            if (params%enable_energy) then
               call update_transport_properties(mesh, flow_mpi, params, transport=transport, T_state=energy%T)
            else
               call update_transport_properties(mesh, flow_mpi, params, transport=transport)
            end if
         end if
         call profiler_stop('Transport_Update')
      end if
      
      ! C. Advance Momentum & Pressure (Projection Method).
      call profiler_start('Projection_Step')
      if (params%enable_variable_density .and. allocated(fields%projection_rho)) then
         fields%projection_rho = transport%rho
      end if
      call advance_projection_step(mesh, flow_mpi, bc, params, transport, fields, stats)
      call profiler_stop('Projection_Step')
      
      ! D. Advance Species Transport.
      if (params%enable_species) then
         call profiler_start('Species_Transport')
         call advance_species_transport(mesh, flow_mpi, bc, params, fields, species, transport)
         call profiler_stop('Species_Transport')
      end if
      
      ! Unified Physics Synchronization Step for load-balancing
      if (params%enable_chemistry_load_balancing) then
         is_chem_step = params%enable_reactions .and. &
                        (mod(step - start_step, params%chemistry_update_interval) == 0 .or. step == params%nsteps)
         is_rad_step = params%enable_energy .and. params%enable_radiation .and. &
                       (mod(step - 1, params%radiation_update_interval) == 0 .or. step == 1)
         if (is_chem_step .or. is_rad_step) then
            call profiler_start('Unified_Physics_Sync')
            call flow_allgather_owned_scalar(flow_mpi, energy%T, energy%T)
            call flow_allgather_owned_matrix(flow_mpi, species%Y, species%Y)
            call profiler_stop('Unified_Physics_Sync')
         end if
      end if
      
      ! E. Optional cell-local Cantera chemistry after species transport.
      ! When chemistry_update_interval > 1, integrate over the accumulated
      ! skipped flow time instead of only one CFD dt.
      if (params%enable_reactions) then
         chemistry_accum_dt = chemistry_accum_dt + params%dt
         if (mod(step - start_step, params%chemistry_update_interval) == 0 .or. step == params%nsteps) then
            call profiler_start('Chemistry_Update')
            call advance_cantera_chemistry(mesh, flow_mpi, params, species%Y, energy, step, &
                                           chemistry_dt=chemistry_accum_dt, &
                                           max_dT_chem=chem_max_dT, max_dY_chem=chem_max_dY, &
                                           max_rel_drho_chem=chem_max_rel_drho, min_source_alpha=chem_min_alpha, &
                                           raw_max_dT_chem=chem_raw_max_dT, raw_max_dY_chem=chem_raw_max_dY, &
                                           raw_max_rel_drho_chem=chem_raw_max_rel_drho)
            chemistry_accum_dt = zero
            call profiler_stop('Chemistry_Update')
         end if
      end if

      ! F. Optional spectral radiation update before energy transport.  On
      ! non-update steps this call preserves the existing qrad field so the
      ! last radiation source is reused.
      if (params%enable_energy .and. params%enable_radiation) then
         if (params%enable_species) then
            call radiation_update_source(mesh, bc, flow_mpi, rad_mpi, rad_context, params, energy, species%Y, step, time, params%dt)
         else
            call radiation_update_source(mesh, bc, flow_mpi, rad_mpi, rad_context, params, energy, step=step, time=time, dt=params%dt)
         end if
      end if

      ! G. Advance passive/reacting sensible-enthalpy transport.
      if (params%enable_energy) then
         call profiler_start('Energy_Transport')
         if (params%enable_variable_density) then
            call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo)
         end if
         if (params%enable_species) then
            call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy, species%Y, transport, step)
         else
            call advance_energy_transport(mesh, flow_mpi, bc, params, fields, energy, transport=transport, step=step)
         end if

         ! No-op for enable_variable_density=.false.; guarded variable-density
         ! mode refreshes transport%rho after the thermo sync.
         call sync_active_density_from_thermo(mesh, flow_mpi, params, transport, energy%rho_thermo)
         if (params%enable_species) then
            call update_lowmach_divergence_source_from_density(mesh, flow_mpi, bc, params, transport, fields, species%Y, energy%T)
         else
            call update_lowmach_divergence_source_from_density(mesh, flow_mpi, bc, params, transport, fields, T_state=energy%T)
         end if

         ! Post-density corrector projection (no-op unless enable_density_corrector_projection=.true.).
         if (params%enable_density_corrector_projection) then
            call profiler_start('Density_Corrector_Projection')
            call correct_projection_to_lowmach_source(mesh, flow_mpi, bc, params, transport, fields, stats)
            call profiler_stop('Density_Corrector_Projection')
         end if

         call profiler_stop('Energy_Transport')
      end if
      
      call compute_step_health(mesh, flow_mpi, params, fields, transport, step_ke_after, step_cfl_after, step_umax_after)
      stats%cfl = step_cfl_after
      stats%kinetic_energy = step_ke_after
      stats%max_velocity = step_umax_after
      call evaluate_step_rejection(params, stats, chem_max_dT, chem_max_dY, chem_max_rel_drho, &
                                   step_ke_before, step_ke_after, step_cfl_after, fields, energy, transport, &
                                   reject_step, reject_reason)
      if (reject_step) then
         step_trial_dt = params%dt
         step_retry_count = step_retry_count + 1
         call restore_step_backup(step_backup, fields, species, energy, transport, stats, time, params%dt, chemistry_accum_dt)
         fields%rhs_old_valid = .false.
         if (params%enable_species) species%rhs_old_valid = .false.
         if (params%enable_energy) energy%rhs_old_valid = .false.
         step_cut_dt = max(max(params%step_reject_min_dt, params%min_dt), step_trial_dt * params%step_reject_cut_factor)
         params%dt = min(step_cut_dt, params%max_dt)
         if (flow_mpi%rank == 0) then
            write(output_unit,'(a,i0,a,i0,a,es12.5,a,es12.5,a,a)') &
               'step rejection: step=', step, ' retry=', step_retry_count, &
               ' old_dt=', step_trial_dt, ' new_dt=', params%dt, ' reason=', trim(reject_reason)
         end if
         if (step_retry_count > params%max_step_retries) then
            call fatal_error('main', 'maximum step retries exceeded at step '//trim(int_to_string(step))//': '//trim(reject_reason))
         end if
         if (params%stop_on_min_dt_failure .and. params%dt <= params%step_reject_min_dt * (1.0_rk + 10.0_rk*tiny_safe)) then
            call fatal_error('main', 'step rejected at minimum dt at step '//trim(int_to_string(step))//': '//trim(reject_reason))
         end if
         cycle
      end if
      step_accepted = .true.
      end do

      time = time + params%dt
      params%dt_old = params%dt
      if (params%enable_reactions) then
         call apply_chemistry_timestep_control(flow_mpi, params, chem_max_dT, chem_max_dY, &
                                               chem_max_rel_drho, chem_min_alpha)
      end if

      call write_restart_file(mesh, flow_mpi, params, fields, species, energy, transport, step, time, chemistry_accum_dt)

      ! E. Periodic terminal status, diagnostics, and file output.
      if (do_status_step) then
         stats%wall_time = real(mpi_wtime(), rk) - t0
         call profiler_start('Flow_Diagnostics')
         call compute_global_observables(mesh, flow_mpi, fields, species, transport, stats)
         call profiler_stop('Flow_Diagnostics')

         if (do_output_step) then
            call profiler_start('Diagnostics_Write_Flow')
            call write_diagnostics_row(params, flow_mpi, step, time, stats)
            call write_boundary_mass_balance_diagnostics(mesh, flow_mpi, params, fields, transport, step, time)
            if (params%enable_variable_density .and. step > 0 .and. time > 0.0_rk) then
               call write_variable_density_diagnostics(mesh, flow_mpi, params, fields, energy, transport, step, time)
            end if
            call profiler_stop('Diagnostics_Write_Flow')

            call profiler_start('Output_Write_VTU')
            call write_vtu_unstructured(params, flow_mpi, mesh, fields, species, energy, transport, step, time)
            call profiler_stop('Output_Write_VTU')
         end if

         if (params%enable_energy .and. do_output_step) then
            call profiler_start('Diagnostics_Write_Energy')
            call write_energy_diagnostics_row(mesh, flow_mpi, params, energy, step, time, &
                                              write_file=do_output_step)
            call profiler_stop('Diagnostics_Write_Energy')
         end if
      end if

      if (do_terminal_step) then
            stats%wall_time = real(mpi_wtime(), rk) - t0

            select case (trim(params%terminal_detail))
            case ('brief')
               if (flow_mpi%rank == 0) then
                  write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,i0,2x,a,a,2x,a,a,2x,a,a)') &
                     'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, &
                     'piter=', stats%pressure_iterations, 'out=', yn(do_output_step), &
                     'rst=', yn(is_restart_write_step(params, step)), &
                     'rad=', yn(is_radiation_update_step(params, step))
               end if

            case ('health')
               call compute_terminal_health(mesh, flow_mpi, params, fields, transport, energy, &
                                            terminal_projection_metric, terminal_ke, terminal_tmin, terminal_tmax)

               if (flow_mpi%rank == 0) then
                  if (params%enable_variable_density) then
                     if (params%enable_energy) then
                        write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') &
                           'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, &
                           'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, &
                           'lm_res_proj=', terminal_projection_metric, 'Tmin=', terminal_tmin, &
                           'Tmax=', terminal_tmax, 'KE=', terminal_ke, &
                           'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), &
                           'rad=', yn(is_radiation_update_step(params, step))
                     else
                        write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') &
                           'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, &
                           'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, &
                           'lm_res_proj=', terminal_projection_metric, 'KE=', terminal_ke, &
                           'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), &
                           'rad=', yn(is_radiation_update_step(params, step))
                     end if
                  else
                     if (params%enable_energy) then
                        write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') &
                           'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, &
                           'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, &
                           'max_div=', terminal_projection_metric, 'Tmin=', terminal_tmin, &
                           'Tmax=', terminal_tmax, 'KE=', terminal_ke, &
                           'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), &
                           'rad=', yn(is_radiation_update_step(params, step))
                     else
                        write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es10.3,2x,a,i0,2x,a,es12.5,2x,a,es12.5,2x,a,a,2x,a,a,2x,a,a)') &
                           'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, &
                           'cfl=', stats%cfl, 'piter=', stats%pressure_iterations, &
                           'max_div=', terminal_projection_metric, 'KE=', terminal_ke, &
                           'out=', yn(do_output_step), 'rst=', yn(is_restart_write_step(params, step)), &
                           'rad=', yn(is_radiation_update_step(params, step))
                     end if
                  end if
               end if

            case ('full')
               if (params%enable_variable_density) then
                  call compute_variable_density_console_diagnostics(mesh, flow_mpi, params, fields, transport, time, &
                                                                    console_max_divu, console_sproj_max, &
                                                                    console_lm_res_projection, console_lm_res_current, &
                                                                    console_mass_balance_defect)
               end if

               if (flow_mpi%rank == 0) then
                  if (params%enable_variable_density) then
                     write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,i0,2x,a,es12.5)') &
                        'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, 'max_divu=', console_max_divu, &
                        'lm_res_proj=', console_lm_res_projection, 'lm_res_cur=', console_lm_res_current, &
                        'Sproj_max=', console_sproj_max, 'mass_def=', console_mass_balance_defect, &
                        'piter=', stats%pressure_iterations, '|U|max=', stats%max_velocity
                  else
                     write(output_unit,'(a,i0,2x,a,es12.5,2x,a,f8.2,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,es12.5,2x,a,i0,2x,a,es12.5)') &
                        'step=', step, 'time=', time, 'clock=', stats%wall_time, 'dt=', params%dt, 'max_div=', stats%max_divergence, &
                        'cfl=', stats%cfl, 'ke=', stats%kinetic_energy, 'piter=', stats%pressure_iterations, '|U|max=', stats%max_velocity
                  end if
               end if
            end select
      end if
   end do

   ! 8. Clean up and finalize simulation.
   call finalize_bc_set(bc)
   call finalize_fields(fields)
   if (params%enable_species) then
      call finalize_species(species)
   end if
   if (params%enable_energy) then
      call finalize_energy(energy)
   end if
   call finalize_transport(transport)
   call mesh_finalize(mesh)

   call profiler_stop('Total_Simulation')
   call profiler_report(flow_mpi%comm, flow_mpi%rank, flow_mpi%nprocs)
   if (params%write_cantera_cache_stats) then
      call write_chemistry_screening_stats(flow_mpi, params)
      call write_cantera_cache_stats(flow_mpi)
   end if
   call radiation_finalize(rad_context)
   call radiation_mpi_finalize(rad_mpi)
   call flow_mpi_finalize(flow_mpi)
   call mpi_flow_shutdown()

contains

   subroutine save_step_backup(fields, species, energy, transport, stats, time, dt, chemistry_accum_dt, backup)
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(in) :: species
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      type(solver_stats_t), intent(in) :: stats
      real(rk), intent(in) :: time, dt, chemistry_accum_dt
      type(step_backup_t), intent(inout) :: backup

      backup%fields = fields
      backup%species = species
      backup%energy = energy
      backup%transport = transport
      backup%stats = stats
      backup%time = time
      backup%dt = dt
      backup%chemistry_accum_dt = chemistry_accum_dt
      backup%valid = .true.
   end subroutine save_step_backup


   subroutine restore_step_backup(backup, fields, species, energy, transport, stats, time, dt, chemistry_accum_dt)
      type(step_backup_t), intent(in) :: backup
      type(flow_fields_t), intent(inout) :: fields
      type(species_fields_t), intent(inout) :: species
      type(energy_fields_t), intent(inout) :: energy
      type(transport_properties_t), intent(inout) :: transport
      type(solver_stats_t), intent(inout) :: stats
      real(rk), intent(inout) :: time, dt, chemistry_accum_dt

      if (.not. backup%valid) call fatal_error('main', 'attempted to restore an invalid step backup')
      fields = backup%fields
      species = backup%species
      energy = backup%energy
      transport = backup%transport
      stats = backup%stats
      time = backup%time
      dt = backup%dt
      chemistry_accum_dt = backup%chemistry_accum_dt
   end subroutine restore_step_backup


   subroutine compute_step_health(mesh, flow, params, fields, transport, kinetic_energy, cfl_value, max_velocity)
      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
      real(rk), intent(out) :: kinetic_energy
      real(rk), intent(out) :: cfl_value
      real(rk), intent(out) :: max_velocity

      integer :: c, lf, f, ierr
      real(rk) :: rho_cell, vmag, local_ke, global_ke
      real(rk) :: local_umax, global_umax
      real(rk) :: local_cfl_rate, global_cfl_rate, cell_outward_flux

      local_ke = zero
      local_umax = zero
      local_cfl_rate = zero

      do c = flow%first_cell, flow%last_cell
         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_safe)
         end if
         vmag = sqrt(dot_product(fields%u(:, c), fields%u(:, c)))
         local_umax = max(local_umax, vmag)
         local_ke = local_ke + 0.5_rk * rho_cell * mesh%cells(c)%volume * vmag * vmag

         if (allocated(fields%face_flux)) then
            cell_outward_flux = zero
            do lf = 1, mesh%ncell_faces(c)
               f = mesh%cell_faces(lf, c)
               cell_outward_flux = cell_outward_flux + abs(fields%face_flux(f))
            end do
            if (mesh%cells(c)%volume > tiny_safe) then
               local_cfl_rate = max(local_cfl_rate, 0.5_rk * cell_outward_flux / mesh%cells(c)%volume)
            end if
         end if
      end do

      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 step kinetic energy')
      call MPI_Allreduce(local_umax, global_umax, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing step max velocity')
      call MPI_Allreduce(local_cfl_rate, global_cfl_rate, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing step CFL rate')

      kinetic_energy = global_ke
      max_velocity = global_umax
      cfl_value = global_cfl_rate * params%dt
   end subroutine compute_step_health


   !> Evaluate multi-guard physical and numerical safety gates for step rejection and retry.
   !!
   !! Evaluates multiple safety constraints over the trial timestep. If any limits are
   !! violated (e.g. NaN/Inf presence, chemistry increments, density vacuum floors,
   !! temperature excursions, or PCG linear solver convergence bounds), the trial step
   !! is rejected, triggering a temporal step reduction and retry.
   !!
   !! @param params Input case parameters and user namelist settings.
   !! @param stats Cumulative solver and PCG convergence metrics.
   !! @param max_dT_chem Maximum temperature increment during the chemistry advance [K].
   !! @param max_dY_chem Maximum species mass fraction increment during the chemistry advance.
   !! @param max_rel_drho_chem Maximum relative density change during the chemistry advance.
   !! @param ke_before Kinetic energy of the domain before the trial timestep [J].
   !! @param ke_after Trial domain kinetic energy after the trial timestep [J].
   !! @param cfl_after Peak domain CFL value computed at the trial state.
   !! @param fields Struct containing current flow fields (u, p).
   !! @param energy Struct containing current thermodynamic/energy fields (T, h).
   !! @param transport Struct containing current transport/thermo properties (rho, mu).
   !! @param reject Output logical indicating if the trial step must be rejected.
   !! @param reason Output character string describing the first tripped safety gate.
   subroutine evaluate_step_rejection(params, stats, max_dT_chem, max_dY_chem, max_rel_drho_chem, &
                                      ke_before, ke_after, cfl_after, fields, energy, transport, &
                                      reject, reason)
      use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
      type(case_params_t), intent(in) :: params
      type(solver_stats_t), intent(in) :: stats
      real(rk), intent(in) :: max_dT_chem, max_dY_chem, max_rel_drho_chem
      real(rk), intent(in) :: ke_before, ke_after, cfl_after
      type(flow_fields_t), intent(in) :: fields
      type(energy_fields_t), intent(in) :: energy
      type(transport_properties_t), intent(in) :: transport
      logical, intent(out) :: reject
      character(len=*), intent(out) :: reason

      real(rk) :: ke_growth, cfl_limit

      reject = .false.
      reason = 'accepted'
      if (.not. params%enable_step_rejection) return

      ! Option B.1: NaN and Infinity Sanitation
      if (params%reject_on_nan) then
         if (any(ieee_is_nan(fields%u)) .or. any(abs(fields%u) > 1.0e8_rk)) then
            reject = .true.
            reason = 'NaN or Inf detected in velocity field u'
            return
         end if
         if (any(ieee_is_nan(fields%p)) .or. any(abs(fields%p) > 1.0e12_rk)) then
            reject = .true.
            reason = 'NaN or Inf detected in pressure field p'
            return
         end if
         if (allocated(energy%T)) then
            if (any(ieee_is_nan(energy%T)) .or. any(energy%T > 1.0e5_rk)) then
               reject = .true.
               reason = 'NaN or Inf detected in temperature field T'
               return
            end if
         end if
         if (allocated(transport%rho)) then
            if (any(ieee_is_nan(transport%rho))) then
               reject = .true.
               reason = 'NaN detected in density field rho'
               return
            end if
         end if
      end if

      ! Option B.2: Chemistry Limits
      if (params%chemistry_max_dT_per_step > zero .and. max_dT_chem > params%chemistry_max_dT_per_step) then
         reject = .true.
         reason = 'chemistry dT exceeded committed limit'
         return
      end if
      if (params%chemistry_max_dY_per_step > zero .and. max_dY_chem > params%chemistry_max_dY_per_step) then
         reject = .true.
         reason = 'chemistry dY exceeded committed limit'
         return
      end if
      if (params%chemistry_max_rel_rho_change_per_step > zero .and. &
          max_rel_drho_chem > params%chemistry_max_rel_rho_change_per_step) then
         reject = .true.
         reason = 'chemistry density change exceeded committed limit'
         return
      end if

      ! Option B.3: CFL Overshoot
      if (params%max_cfl_overshoot_factor > zero) then
         cfl_limit = params%max_cfl * params%max_cfl_overshoot_factor
         if (cfl_after > cfl_limit) then
            reject = .true.
            reason = 'post-step CFL overshoot'
            return
         end if
      end if

      ! Option B.4: Kinetic Energy Growth
      if (params%max_KE_growth_per_step > zero .and. ke_before > params%ke_reject_floor) then
         ke_growth = ke_after / ke_before
         if (ke_growth > params%max_KE_growth_per_step) then
            reject = .true.
            reason = 'kinetic energy growth exceeded limit'
            return
         end if
      end if

      ! Option B.5: PCG Pressure Iterations Ceiling
      if (params%reject_on_pressure_max_iter .and. stats%pressure_iterations >= params%pressure_max_iter) then
         reject = .true.
         reason = 'pressure solve reached max iterations'
         return
      end if

      ! Option B.6: PCG Convergence Health Audit
      if (stats%pressure_residual > params%max_pressure_residual_limit .and. &
          stats%pressure_residual_abs > params%pressure_abs_tol) then
         reject = .true.
         reason = 'pressure solver residual exceeds stability tolerance'
         return
      end if

      ! Option B.7: Density Floor/Vacuum Guard
      if (allocated(transport%rho)) then
         if (minval(transport%rho) < params%rho * params%min_density_ratio_limit) then
            reject = .true.
            reason = 'density vacuum limit violated'
            return
         end if
      end if

      ! Option B.8: Temperature Floor and Ceiling Excursions
      if (allocated(energy%T)) then
         if (minval(energy%T) < params%min_temperature_limit) then
            reject = .true.
            reason = 'temperature fell below safety floor'
            return
         end if
         if (maxval(energy%T) > params%max_temperature_limit) then
            reject = .true.
            reason = 'temperature exceeded safety ceiling'
            return
         end if
      end if
   end subroutine evaluate_step_rejection


   function int_to_string(value) result(text)
      integer, intent(in) :: value
      character(len=32) :: text
      write(text,'(i0)') value
   end function int_to_string


   !> Cut the next dynamic timestep when chemistry produced a large state jump.
   !! This does not reject the current step; it prevents the next operator-split
   !! chemistry/projection update from repeating an overly large heat-release jump.
   subroutine apply_chemistry_timestep_control(flow, params, max_dT, max_dY, max_rel_drho, min_alpha)
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(inout) :: params
      real(rk), intent(in) :: max_dT, max_dY, max_rel_drho, min_alpha

      real(rk) :: factor, old_dt, threshold

      if (.not. params%use_dynamic_dt) return
      factor = one

      threshold = params%chemistry_max_dT_per_step
      if (threshold > zero .and. max_dT > threshold) then
         factor = min(factor, params%chemistry_dt_safety * threshold / max(max_dT, tiny(1.0_rk)))
      end if

      threshold = params%chemistry_max_dY_per_step
      if (threshold > zero .and. max_dY > threshold) then
         factor = min(factor, params%chemistry_dt_safety * threshold / max(max_dY, tiny(1.0_rk)))
      end if

      threshold = params%chemistry_max_rel_rho_change_per_step
      if (threshold > zero .and. max_rel_drho > threshold) then
         factor = min(factor, params%chemistry_dt_safety * threshold / max(max_rel_drho, tiny(1.0_rk)))
      end if

      factor = max(params%chemistry_dt_min_factor, min(one, factor))
      if (factor < one) then
         old_dt = params%dt
         params%dt = max(params%min_dt, min(params%max_dt, params%dt * factor))
         if (flow%rank == 0) then
            write(output_unit,'(a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4,a,es11.4)') &
               'chemistry dt control: dt ', old_dt, ' -> ', params%dt, &
               ' max_dT=', max_dT, ' max_dY=', max_dY, &
               ' max_rel_drho=', max_rel_drho, ' min_alpha=', min_alpha
         end if
      end if
   end subroutine apply_chemistry_timestep_control

   !> Apply a one-time spherical ignition kernel during initialization.
   !!
   !! The current first version supports one spherical overwrite kernel.  It is
   !! skipped automatically for restart runs so restart state is not modified by
   !! fresh-run initialization controls.
   subroutine apply_ignition_kernel(mesh, flow, params, species, energy)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(case_params_t), intent(in) :: params
      type(species_fields_t), intent(inout) :: species
      type(energy_fields_t), intent(inout) :: energy

      integer :: c, affected_cells
      real(rk) :: delta(3), radius2
      real(rk) :: kernel_Y(max_species)
      logical :: set_composition

      if (.not. params%enable_ignition_kernel) return
      if (.not. params%enable_energy) return
      if (.not. energy%initialized) return

      affected_cells = 0
      radius2 = params%ignition_kernel_radius * params%ignition_kernel_radius
      set_composition = .false.
      kernel_Y = zero

      if (params%enable_species .and. species%nspecies > 0 .and. &
          composition_string_is_set(params%ignition_kernel_composition)) then
         call parse_named_composition(params%ignition_kernel_composition, species%names, species%nspecies, &
                                      kernel_Y, 'ignition_kernel_composition', &
                                      normalize=.true., require_sum_positive=.true.)
         set_composition = .true.
      end if

      do c = 1, mesh%ncells
         delta = mesh%cells(c)%center - params%ignition_kernel_center
         if (dot_product(delta, delta) <= radius2) then
            energy%T(c) = params%ignition_kernel_T
            if (set_composition) species%Y(1:species%nspecies, c) = kernel_Y(1:species%nspecies)
            affected_cells = affected_cells + 1
         end if
      end do

      if (set_composition) then
         species%Y_old = species%Y
         species%rhs_old = zero
         species%rhs_old_valid = .false.
      end if

      if (flow%rank == 0) then
         write(output_unit,'(a,i0,a,es12.5,a,3(1x,es12.5))') &
            'Ignition kernel initialized: cells=', affected_cells, &
            ' T=', params%ignition_kernel_T, ' center=', params%ignition_kernel_center
      end if
   end subroutine apply_ignition_kernel



   character(len=1) function yn(flag) result(mark)
      logical, intent(in) :: flag
      if (flag) then
         mark = 'Y'
      else
         mark = 'N'
      end if
   end function yn

   logical function is_restart_write_step(params, step) result(write_now)
      type(case_params_t), intent(in) :: params
      integer, intent(in) :: step

      write_now = .false.
      if (.not. params%write_restart) return
      if (params%restart_interval <= 0) return
      write_now = (mod(step, params%restart_interval) == 0) .or. (step == params%nsteps)
   end function is_restart_write_step

   logical function is_radiation_update_step(params, step) result(update_now)
      type(case_params_t), intent(in) :: params
      integer, intent(in) :: step

      update_now = .false.
      if (.not. params%enable_energy) return
      if (.not. params%enable_radiation) return
      if (params%radiation_update_interval <= 0) return
      update_now = (mod(step - 1, params%radiation_update_interval) == 0) .or. (step == 1)
   end function is_radiation_update_step

   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

   !> Write a concise, reproducible summary of physics and numerical settings.
   subroutine write_run_config_summary(mesh, flow, rad, params, case_file)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(radiation_mpi_t), intent(in) :: rad
      type(case_params_t), intent(in) :: params
      character(len=*), intent(in) :: case_file

      integer :: unit_id, ios
      character(len=1024) :: filename

      if (flow%rank /= 0) return

      call write_run_config_summary_to_unit(output_unit, mesh, flow, rad, params, case_file)

      filename = trim(params%output_dir) // '/run_config_summary.txt'
      open(newunit=unit_id, file=trim(filename), status='replace', action='write', iostat=ios)
      if (ios /= 0) then
         write(output_unit,'(a,a)') 'warning: could not write run configuration summary: ', trim(filename)
         return
      end if
      call write_run_config_summary_to_unit(unit_id, mesh, flow, rad, params, case_file)
      close(unit_id)
   end subroutine write_run_config_summary


   subroutine write_run_config_summary_to_unit(unit_id, mesh, flow, rad, params, case_file)
      integer, intent(in) :: unit_id
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(radiation_mpi_t), intent(in) :: rad
      type(case_params_t), intent(in) :: params
      character(len=*), intent(in) :: case_file

      integer :: i

      write(unit_id,'(a)') '======================================================================'
      write(unit_id,'(a)') 'RUN CONFIGURATION SUMMARY'
      write(unit_id,'(a)') '======================================================================'
      write(unit_id,'(a,a)') 'case: ', trim(case_file)
      write(unit_id,'(a,a)') 'mesh_dir: ', trim(params%mesh_dir)
      write(unit_id,'(a,i0)') 'cells: ', mesh%ncells
      write(unit_id,'(a,i0)') 'faces: ', mesh%nfaces
      write(unit_id,'(a,i0)') 'flow/chemistry MPI ranks: ', flow%nprocs
      write(unit_id,'(a,i0)') 'radiation MPI ranks: ', rad%nprocs
      write(unit_id,'(a,i0,a,es12.5,a,i0,a,i0,a,l1)') 'time: nsteps=', params%nsteps, ' dt=', params%dt, &
         ' output_interval=', params%output_interval, ' terminal_interval=', params%terminal_interval, &
         ' terminal_detail='//trim(params%terminal_detail)//' dynamic_dt=', params%use_dynamic_dt
      write(unit_id,'(a,es12.5,a,es12.5,a,es12.5)') 'time limits: max_dt=', params%max_dt, &
         ' min_dt=', params%min_dt, ' dt_growth_limit=', params%dt_growth_limit
      write(unit_id,'(a,l1)')      '  step rejection enabled           : ', params%enable_step_rejection
      write(unit_id,'(a,i0)')      '  max step retries                 : ', params%max_step_retries
      write(unit_id,'(a,es12.5)') '  step reject cut factor           : ', params%step_reject_cut_factor
      write(unit_id,'(a,es12.5)') '  step reject minimum dt           : ', params%step_reject_min_dt
      write(unit_id,'(a,es12.5)') '  max KE growth per step           : ', params%max_KE_growth_per_step
      write(unit_id,'(a,es12.5)') '  KE reject floor                  : ', params%ke_reject_floor
      write(unit_id,'(a,es12.5)') '  max CFL overshoot factor         : ', params%max_cfl_overshoot_factor
      write(unit_id,'(a,l1)')      '  reject on pressure max iter      : ', params%reject_on_pressure_max_iter
      write(unit_id,'(a,l1)')      '  stop on min-dt failure           : ', params%stop_on_min_dt_failure

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Numerical schemes:'
      write(unit_id,'(a,a)') '  legacy convection_scheme     : ', trim(params%convection_scheme)
      write(unit_id,'(a,a)') '  momentum convection          : ', trim(params%momentum_convection_scheme)
      write(unit_id,'(a,a)') '  species convection           : ', trim(params%species_convection_scheme)
      write(unit_id,'(a,a)') '  energy convection            : ', trim(params%energy_convection_scheme)
      write(unit_id,'(a,a)') '  scalar limiter               : ', trim(params%scalar_limiter)
      write(unit_id,'(a,a)') '  species time scheme          : ', trim(params%species_time_scheme)
      write(unit_id,'(a,a)') '  energy time scheme           : ', trim(params%energy_time_scheme)
      write(unit_id,'(a,i0)') '  pressure max iterations      : ', params%pressure_max_iter
      write(unit_id,'(a,es12.5)') '  pressure tolerance           : ', params%pressure_tol

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Flow and thermodynamics:'
      write(unit_id,'(a,l1)') '  variable density             : ', params%enable_variable_density
      write(unit_id,'(a,a)') '  density EOS                  : ', trim(params%density_eos)
      write(unit_id,'(a,l1)') '  variable viscosity           : ', params%enable_variable_nu
      write(unit_id,'(a,i0)') '  transport update interval    : ', params%transport_update_interval
      write(unit_id,'(a,l1)') '  energy equation              : ', params%enable_energy
      write(unit_id,'(a,l1)') '  Cantera thermo               : ', params%enable_cantera_thermo
      write(unit_id,'(a,l1)') '  species enthalpy diffusion   : ', params%enable_species_enthalpy_diffusion
      write(unit_id,'(a,l1)') '  upwind dilatation density    : ', params%upwind_dilatation_density
      write(unit_id,'(a,es12.5)') '  background pressure [Pa]     : ', params%background_press
      write(unit_id,'(a,es12.5)') '  background temperature [K]   : ', params%background_temp

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Chemistry timestep/source controls:'
      write(unit_id,'(a,es12.5)') '  chemistry max dT per step [K]      : ', params%chemistry_max_dT_per_step
      write(unit_id,'(a,es12.5)') '  chemistry max rel rho change       : ', params%chemistry_max_rel_rho_change_per_step
      write(unit_id,'(a,es12.5)') '  chemistry max dY per step          : ', params%chemistry_max_dY_per_step
      write(unit_id,'(a,es12.5)') '  chemistry dt safety                : ', params%chemistry_dt_safety
      write(unit_id,'(a,l1)')      '  chemistry source limiter           : ', params%chemistry_limit_source_update
      write(unit_id,'(a,es12.5)') '  chemistry source relaxation        : ', params%chemistry_source_relaxation

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Species and chemistry:'
      write(unit_id,'(a,l1)') '  species transport            : ', params%enable_species
      write(unit_id,'(a,i0)') '  transported species count    : ', params%nspecies
      write(unit_id,'(a,l1)') '  reactions enabled            : ', params%enable_reactions
      write(unit_id,'(a,a)') '  Cantera mechanism            : ', trim(params%cantera_mech_file)
      if (len_trim(params%cantera_phase_name) > 0) then
         write(unit_id,'(a,a)') '  Cantera phase                : ', trim(params%cantera_phase_name)
      else
         write(unit_id,'(a)') '  Cantera phase                : <default first phase>'
      end if
      write(unit_id,'(a,i0)') '  Cantera reaction count       : ', params%cantera_reaction_count
      write(unit_id,'(a,i0)') '  chemistry update interval    : ', params%chemistry_update_interval
      write(unit_id,'(a,es12.5)') '  chemistry rtol               : ', params%chemistry_rtol
      write(unit_id,'(a,es12.5)') '  chemistry atol               : ', params%chemistry_atol
      write(unit_id,'(a,i0)') '  chemistry max steps          : ', params%chemistry_max_steps
      write(unit_id,'(a,l1)') '  chemistry energy enabled     : ', params%chemistry_energy_enabled
      write(unit_id,'(a,es12.5)') '  chemistry T cutoff [K]       : ', params%chemistry_temperature_cutoff
      write(unit_id,'(a,es12.5)') '  chemistry min reactive Y     : ', params%chemistry_min_reactive_mass_fraction
      write(unit_id,'(a,l1)') '  chemistry_active_screening   : ', &
         params%chemistry_active_species_threshold > zero .and. params%chemistry_n_active_species > 0
      write(unit_id,'(a,i0)') '  chemistry_n_active_species   : ', params%chemistry_n_active_species
      write(unit_id,'(a,es12.5)') '  chemistry_active_species_threshold : ', params%chemistry_active_species_threshold
      if (params%chemistry_n_active_species > 0) then
         do i = 1, params%chemistry_n_active_species
            if (len_trim(params%chemistry_active_species_name(i)) > 0) then
               write(unit_id,'(a,i0,a,a)') '  chemistry_active_species_name(', i, ') : ', &
                  trim(params%chemistry_active_species_name(i))
            end if
         end do
      end if
      if (len_trim(params%initial_composition) > 0) then
         write(unit_id,'(a,a)') '  initial composition          : ', trim(params%initial_composition)
      end if

      call write_ignition_config_summary_to_unit(unit_id, mesh, params)

      call write_boundary_config_summary_to_unit(unit_id, params)

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Radiation coupling:'
      write(unit_id,'(a,l1)') '  radiation enabled            : ', params%enable_radiation
      write(unit_id,'(a,a)') '  radiation pressure source    : ', trim(params%radiation_pressure_source)
      write(unit_id,'(a,a)') '  radiation model              : ', trim(params%radiation_source_model)
      write(unit_id,'(a,i0)') '  radiation update interval    : ', params%radiation_update_interval
      write(unit_id,'(a,i0)') '  n wavenumbers                : ', params%radiation_n_wavenumbers
      write(unit_id,'(a,i0,a,i0)') '  local wavenumber range       : ', rad%first_wavenumber, '-', rad%last_wavenumber
      write(unit_id,'(a,i0)') '  selected radiation species   : ', params%radiation_n_species
      write(unit_id,'(a,l1)') '  radiation debug checks       : ', params%radiation_debug
      write(unit_id,'(a,l1)') '  radiation diagnostics        : ', params%write_radiation_diagnostics

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Restart:'
      write(unit_id,'(a,l1)') '  restart from file            : ', params%restart_from_file
      write(unit_id,'(a,a)') '  restart file                 : ', trim(params%restart_file)
      write(unit_id,'(a,l1)') '  write restart                : ', params%write_restart
      write(unit_id,'(a,i0)') '  restart interval             : ', params%restart_interval
      write(unit_id,'(a,a)') '  restart output dir           : ', trim(params%restart_output_dir)
      write(unit_id,'(a,i0)') '  restart keep                 : ', params%restart_keep
      write(unit_id,'(a,l1)') '  restart use file dt          : ', params%restart_use_file_dt
      write(unit_id,'(a,l1)') '  restart nsteps additional    : ', params%restart_nsteps_are_additional
      write(unit_id,'(a)') '======================================================================'
   end subroutine write_run_config_summary_to_unit


   !> Write ignition-kernel settings to the run configuration summary.
   subroutine write_ignition_config_summary_to_unit(unit_id, mesh, params)
      integer, intent(in) :: unit_id
      type(mesh_t), intent(in) :: mesh
      type(case_params_t), intent(in) :: params

      integer :: c, affected_cells
      real(rk) :: delta(3), radius2

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Ignition initialization:'
      write(unit_id,'(a,l1)') '  enable_ignition_kernel    : ', params%enable_ignition_kernel
      if (.not. params%enable_ignition_kernel) return

      affected_cells = 0
      radius2 = params%ignition_kernel_radius * params%ignition_kernel_radius
      do c = 1, mesh%ncells
         delta = mesh%cells(c)%center - params%ignition_kernel_center
         if (dot_product(delta, delta) <= radius2) affected_cells = affected_cells + 1
      end do

      write(unit_id,'(a,a)') '  ignition_kernel_shape     : ', trim(params%ignition_kernel_shape)
      write(unit_id,'(a,3(1x,es12.5))') '  ignition_kernel_center    : ', params%ignition_kernel_center
      write(unit_id,'(a,es12.5)') '  ignition_kernel_radius    : ', params%ignition_kernel_radius
      write(unit_id,'(a,es12.5)') '  ignition_kernel_T [K]     : ', params%ignition_kernel_T
      write(unit_id,'(a,a)') '  ignition_kernel_blend     : ', trim(params%ignition_kernel_blend)
      write(unit_id,'(a,i0)') '  ignition_kernel_cells     : ', affected_cells
      if (len_trim(params%ignition_kernel_composition) > 0) then
         write(unit_id,'(a,a)') '  ignition_kernel_composition : ', trim(params%ignition_kernel_composition)
      else
         write(unit_id,'(a)') '  ignition_kernel_composition : <unchanged from initial field>'
      end if
      if (params%restart_from_file) then
         write(unit_id,'(a)') '  note                       : restart_from_file=.true.; fresh-run kernel is skipped'
      end if
   end subroutine write_ignition_config_summary_to_unit


   !> Write boundary-condition input settings to the run configuration summary.
   subroutine write_boundary_config_summary_to_unit(unit_id, params)
      integer, intent(in) :: unit_id
      type(case_params_t), intent(in) :: params

      integer :: p

      write(unit_id,'(a)') ''
      write(unit_id,'(a)') 'Boundary conditions:'
      write(unit_id,'(a,i0)') '  number of patches            : ', params%n_patches
      write(unit_id,'(a)') '  mass-flux velocity convention:'
      write(unit_id,'(a)') '    mass_flux/signed_mass_flux : signed rho*dot(U,n_out) [kg/m^2/s]'
      write(unit_id,'(a)') '    inlet_mass_flux            : positive magnitude directed into the domain'
      write(unit_id,'(a)') '    outlet_mass_flux           : positive magnitude directed out of the domain'

      if (params%n_patches <= 0) return

      do p = 1, params%n_patches
         write(unit_id,'(a,i0,a,a)') '  patch[', p, '] ', trim(params%patch_name(p))
         write(unit_id,'(a,a)') '    patch type                 : ', trim(params%patch_type(p))
         write(unit_id,'(a,a)') '    velocity type              : ', trim(params%patch_velocity_type(p))
         write(unit_id,'(a,3(es12.5,1x))') '    velocity [m/s]             : ', &
            params%patch_u(p), params%patch_v(p), params%patch_w(p)
         write(unit_id,'(a,es12.5)') '    mass flux [kg/m^2/s]       : ', params%patch_mass_flux(p)
         write(unit_id,'(a,a)') '    pressure type              : ', trim(params%patch_pressure_type(p))
         write(unit_id,'(a,es12.5)') '    pressure value             : ', params%patch_p(p)
         write(unit_id,'(a,es12.5)') '    pressure gradient          : ', params%patch_dpdn(p)
         write(unit_id,'(a,a)') '    temperature type           : ', trim(params%patch_temperature_type(p))
         write(unit_id,'(a,es12.5)') '    temperature [K]            : ', params%patch_T(p)
         write(unit_id,'(a,a)') '    species type               : ', trim(params%patch_species_type(p))
         if (len_trim(params%patch_composition(p)) > 0) then
            write(unit_id,'(a,a)') '    composition                : ', trim(params%patch_composition(p))
         end if
      end do
   end subroutine write_boundary_config_summary_to_unit


   !> Computes global integral quantities for diagnostic reporting.
   !!
   !! Performs Allreduce operations to find the maximum velocity magnitude, 
   !! total system mass, and minimum species mass fraction.
   subroutine compute_global_observables(mesh, flow, fields, species, transport, stats)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(in) :: flow
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(in) :: species
      type(transport_properties_t), intent(in) :: transport
      type(solver_stats_t), intent(inout) :: stats

      integer :: c
      real(rk) :: local_max_vel, local_min_y, local_total_mass
      real(rk) :: global_max_vel, global_min_y, global_total_mass
      real(rk) :: u, v, w, vel_mag

      local_max_vel = zero
      local_min_y = huge(one)
      local_total_mass = zero

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

         u = fields%u(1, c)
         v = fields%u(2, c)
         w = fields%u(3, c)
         vel_mag = sqrt(u**2 + v**2 + w**2)
         if (vel_mag > local_max_vel) local_max_vel = vel_mag

         local_total_mass = local_total_mass + (transport%rho(c) * mesh%cells(c)%volume)

         if (species%nspecies > 0) then
            local_min_y = min(local_min_y, minval(species%Y(:, c)))
         end if
      end do

      call MPI_Allreduce(local_max_vel, global_max_vel, 1, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm)
      call MPI_Allreduce(local_total_mass, global_total_mass, 1, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm)
      
      if (species%nspecies > 0) then
         call MPI_Allreduce(local_min_y, global_min_y, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm)
      else
         global_min_y = zero
      end if

      stats%max_velocity = global_max_vel
      stats%total_mass = global_total_mass
      stats%min_species_y = global_min_y

   end subroutine compute_global_observables

   !> Compute concise variable-density low-Mach diagnostics for the terminal.
   !!
   !! In variable-density mode raw div(u) is not a projection-error metric
   !! because the low-Mach target is div(u)=S.  This helper reports both raw
   !! div(u) and residuals against the source actually used by projection.
   subroutine compute_variable_density_console_diagnostics(mesh, flow, params, fields, transport, time, &
                                                           max_divu, max_s_projection, &
                                                           max_lm_res_projection, max_lm_res_current, &
                                                           mass_balance_defect)
      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
      real(rk), intent(in) :: time
      real(rk), intent(out) :: max_divu
      real(rk), intent(out) :: max_s_projection
      real(rk), intent(out) :: max_lm_res_projection
      real(rk), intent(out) :: max_lm_res_current
      real(rk), intent(out) :: mass_balance_defect

      integer :: c, lf, fid, f, owner, nb, ierr
      real(rk) :: vol, divu, flux_out, s_current, s_projection
      real(rk) :: local_max_divu, local_max_s_projection
      real(rk) :: local_max_lm_res_projection, local_max_lm_res_current
      real(rk) :: local_mass, global_mass, local_net_mass_flux, net_mass_flux
      real(rk) :: dt_since_previous, mass_rate
      real(rk) :: local_max_vec(4), global_max_vec(4)
      real(rk) :: local_sum_vec(2), global_sum_vec(2)
      logical, save :: initialized = .false.
      real(rk), save :: previous_mass = 0.0_rk
      real(rk), save :: previous_time = 0.0_rk

      max_divu = 0.0_rk
      max_s_projection = 0.0_rk
      max_lm_res_projection = 0.0_rk
      max_lm_res_current = 0.0_rk
      mass_balance_defect = 0.0_rk

      if (.not. params%enable_variable_density) return
      if (.not. allocated(fields%face_flux)) return

      local_max_divu = 0.0_rk
      local_max_s_projection = 0.0_rk
      local_max_lm_res_projection = 0.0_rk
      local_max_lm_res_current = 0.0_rk
      local_mass = 0.0_rk
      local_net_mass_flux = 0.0_rk

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

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

         divu = 0.0_rk
         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

         s_current = 0.0_rk
         if (allocated(fields%divergence_source)) then
            if (size(fields%divergence_source) >= c) s_current = fields%divergence_source(c)
         end if

         s_projection = s_current
         if (allocated(fields%projection_divergence_source)) then
            if (size(fields%projection_divergence_source) >= c) s_projection = fields%projection_divergence_source(c)
         end if

         local_max_divu = max(local_max_divu, abs(divu))
         local_max_s_projection = max(local_max_s_projection, abs(s_projection))
         local_max_lm_res_projection = max(local_max_lm_res_projection, abs(divu - s_projection))
         local_max_lm_res_current = max(local_max_lm_res_current, abs(divu - s_current))

         if (allocated(transport%rho)) then
            if (size(transport%rho) >= c) local_mass = local_mass + transport%rho(c) * vol
         end if
      end do

      if (allocated(fields%mass_flux)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0) then
               if (owner >= 1 .and. owner <= mesh%ncells) then
                  if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + fields%mass_flux(f)
               end if
            end if
         end do
      else if (allocated(fields%face_flux) .and. allocated(transport%rho)) then
         do f = 1, mesh%nfaces
            owner = mesh%faces(f)%owner
            nb = mesh%faces(f)%neighbor
            if (nb <= 0) nb = mesh%faces(f)%periodic_neighbor

            if (nb <= 0) then
               if (owner >= 1 .and. owner <= mesh%ncells) then
                  if (flow%owned(owner)) local_net_mass_flux = local_net_mass_flux + transport%rho(owner) * fields%face_flux(f)
               end if
            end if
         end do
      end if

      local_max_vec = [local_max_divu, local_max_s_projection, &
                       local_max_lm_res_projection, local_max_lm_res_current]
      call MPI_Allreduce(local_max_vec, global_max_vec, 4, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density max diagnostics')
      max_divu = global_max_vec(1)
      max_s_projection = global_max_vec(2)
      max_lm_res_projection = global_max_vec(3)
      max_lm_res_current = global_max_vec(4)

      local_sum_vec = [local_mass, local_net_mass_flux]
      call MPI_Allreduce(local_sum_vec, global_sum_vec, 2, MPI_DOUBLE_PRECISION, MPI_SUM, flow%comm, ierr)
      if (ierr /= MPI_SUCCESS) call fatal_error('main', 'MPI failure reducing variable-density sum diagnostics')
      global_mass = global_sum_vec(1)
      net_mass_flux = global_sum_vec(2)

      if (initialized) then
         dt_since_previous = time - previous_time
         if (dt_since_previous > tiny(1.0_rk)) then
            mass_rate = (global_mass - previous_mass) / dt_since_previous
            mass_balance_defect = mass_rate + net_mass_flux
         else
            mass_balance_defect = 0.0_rk
         end if
      else
         mass_balance_defect = 0.0_rk
         initialized = .true.
      end if

      previous_mass = global_mass
      previous_time = time

   end subroutine compute_variable_density_console_diagnostics




   !> Parses command line arguments to find the case configuration file.
   !!
   !! If no argument is provided, defaults to `case.nml`.
   subroutine get_case_filename(filename)
      character(len=*), intent(out) :: filename

      integer :: argc

      argc = command_argument_count()
      if (argc >= 1) then
         call get_command_argument(1, filename)
      else
         filename = 'case.nml'
      end if

      if (len_trim(filename) == 0) call fatal_error('main', 'empty case filename')
   end subroutine get_case_filename

end program lowmach_react_hex

