advance_cantera_chemistry Subroutine

   subroutine advance_cantera_chemistry(mesh, flow, params, species_Y, energy, step, chemistry_dt, &
                                        max_dT_chem, max_dY_chem, max_rel_drho_chem, min_source_alpha, &
                                        raw_max_dT_chem, raw_max_dY_chem, raw_max_rel_drho_chem)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(case_params_t), intent(in) :: params
      real(rk), intent(inout) :: species_Y(:,:)
      type(energy_fields_t), intent(inout) :: energy
      integer, intent(in) :: step
      real(rk), intent(in), optional :: chemistry_dt
      real(rk), intent(out), optional :: max_dT_chem, max_dY_chem, max_rel_drho_chem, min_source_alpha
      real(rk), intent(out), optional :: raw_max_dT_chem, raw_max_dY_chem, raw_max_rel_drho_chem

      integer :: nloc, nactive, global_active, i, c, k, n_len, ierr, j, n_sub
      integer :: nactive_screen_species, active_screen_species(max_species)
      integer :: skipped_temperature_local, skipped_named_species_local, skipped_legacy_mass_fraction_local
      integer, allocatable :: active_cell(:)
      real(rk) :: sum_Y, dt_local, reactive_indicator, dt_sub
      real(rk) :: dT_abs, dY_abs, rel_drho, alpha, alpha_base
      real(rk) :: max_dT_local, max_dY_local, max_rel_drho_local, min_alpha_local
      real(rk) :: raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local
      real(rk) :: qdot_max_local, qdot_integral_local
      real(rk) :: reduce_in(3), reduce_out(3), raw_reduce_in(3), raw_reduce_out(3), alpha_reduce
      logical :: use_named_species_screen, use_source_limiter
      real(rk) :: t_total_start, t_reactor_start, total_time_local, reactor_time_local
      real(rk), allocatable :: T_local(:), P_local(:), Y_local(:,:)
      real(rk), allocatable :: T_before(:), rho_before(:), Y_before(:,:)
      real(rk), allocatable :: h_local(:), cp_local(:), lambda_local(:), rho_local(:), alpha_local(:), qdot_local(:)
      character(kind=c_char), allocatable :: c_names_flat(:)
      real(rk), allocatable :: prev_chem_dT_dt(:)
      real(rk) :: dT_max_est

      ! Load-balancing packed dynamic buffers
      integer :: nfields, r
      real(rk) :: qdot_val, T_prev, rho_prev, dY_max_val
      integer, allocatable :: pack_cell(:), all_active_cells(:)
      real(rk), allocatable :: pack_data(:,:), all_data_new(:,:)
      integer, allocatable :: active_counts(:), active_displs(:)
      integer, allocatable :: active_matrix_counts(:), active_matrix_displs(:)

      n_sub = 1
      if (present(max_dT_chem)) max_dT_chem = zero
      if (present(max_dY_chem)) max_dY_chem = zero
      if (present(max_rel_drho_chem)) max_rel_drho_chem = zero
      if (present(min_source_alpha)) min_source_alpha = 1.0_rk
      if (present(raw_max_dT_chem)) raw_max_dT_chem = zero
      if (present(raw_max_dY_chem)) raw_max_dY_chem = zero
      if (present(raw_max_rel_drho_chem)) raw_max_rel_drho_chem = zero

      if (.not. params%enable_reactions) return
      if (params%chemistry_update_interval <= 0) return

      t_total_start = real(MPI_Wtime(), rk)
      reactor_time_local = zero
      total_time_local = zero

      if (.not. params%enable_cantera_thermo) then
         call fatal_error('chemistry', 'Cantera chemistry requires enable_cantera_thermo=.true.')
      end if
      if (.not. energy%initialized) then
         call fatal_error('chemistry', 'energy must be initialized before chemistry')
      end if
      if (params%nspecies <= 0) then
         call fatal_error('chemistry', 'Cantera chemistry requires nspecies > 0 after mechanism discovery')
      end if
      if (size(species_Y,1) < params%nspecies .or. size(species_Y,2) < mesh%ncells) then
         call fatal_error('chemistry', 'species_Y has incompatible shape')
      end if

      call resolve_chemistry_active_species(params, nactive_screen_species, active_screen_species)
      use_named_species_screen = params%chemistry_active_species_threshold > zero .and. &
                                 nactive_screen_species > 0

      if (present(chemistry_dt)) then
         dt_local = chemistry_dt
      else if (step <= 1) then
         dt_local = params%dt
      else
         dt_local = params%dt * real(params%chemistry_update_interval, rk)
      end if
      if (dt_local <= zero) return

      if (allocated(energy%chem_dT_dt)) then
         allocate(prev_chem_dT_dt(size(energy%chem_dT_dt)))
         prev_chem_dT_dt = energy%chem_dT_dt
      end if
      if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate = zero
      if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt = zero
      if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt = zero
      if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt = zero
      if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha = one
      if (allocated(energy%chem_active)) energy%chem_active = zero

      if (params%enable_chemistry_load_balancing) then
         nloc = mesh%ncells
         allocate(active_cell(nloc))
         nactive = 0
         skipped_temperature_local = 0
         skipped_named_species_local = 0
         skipped_legacy_mass_fraction_local = 0
         
         global_active = 0
         do c = 1, mesh%ncells
            if (params%chemistry_temperature_cutoff > zero) then
               if (energy%T(c) < params%chemistry_temperature_cutoff) then
                  skipped_temperature_local = skipped_temperature_local + 1
                  cycle
               end if
            end if
            if (use_named_species_screen) then
               reactive_indicator = zero
               do k = 1, nactive_screen_species
                  reactive_indicator = max(reactive_indicator, &
                                           max(zero, species_Y(active_screen_species(k), c)))
               end do
               if (reactive_indicator < params%chemistry_active_species_threshold) then
                  skipped_named_species_local = skipped_named_species_local + 1
                  cycle
               end if
            end if
            if (params%chemistry_min_reactive_mass_fraction > zero) then
               reactive_indicator = maxval(max(zero, species_Y(1:params%nspecies, c)))
               if (reactive_indicator < params%chemistry_min_reactive_mass_fraction) then
                  skipped_legacy_mass_fraction_local = skipped_legacy_mass_fraction_local + 1
                  cycle
               end if
            end if
            
            global_active = global_active + 1
            if (mod(global_active - 1, flow%nprocs) == flow%rank) then
               nactive = nactive + 1
               active_cell(nactive) = c
            end if
         end do
      else
         nloc = max(0, flow%last_cell - flow%first_cell + 1)
         allocate(active_cell(max(1, nloc)))
         nactive = 0
         skipped_temperature_local = 0
         skipped_named_species_local = 0
         skipped_legacy_mass_fraction_local = 0
         if (nloc > 0) then
            do i = 1, nloc
               c = flow%first_cell + i - 1
               if (params%chemistry_temperature_cutoff > zero) then
                  if (energy%T(c) < params%chemistry_temperature_cutoff) then
                     skipped_temperature_local = skipped_temperature_local + 1
                     cycle
                  end if
               end if
               if (use_named_species_screen) then
                  reactive_indicator = zero
                  do k = 1, nactive_screen_species
                     reactive_indicator = max(reactive_indicator, &
                                              max(zero, species_Y(active_screen_species(k), c)))
                  end do
                  if (reactive_indicator < params%chemistry_active_species_threshold) then
                     skipped_named_species_local = skipped_named_species_local + 1
                     cycle
                  end if
               end if
               if (params%chemistry_min_reactive_mass_fraction > zero) then
                  reactive_indicator = maxval(max(zero, species_Y(1:params%nspecies, c)))
                  if (reactive_indicator < params%chemistry_min_reactive_mass_fraction) then
                     skipped_legacy_mass_fraction_local = skipped_legacy_mass_fraction_local + 1
                     cycle
                  end if
               end if
               nactive = nactive + 1
               active_cell(nactive) = c
            end do
         end if
         
         call MPI_Allreduce(nactive, global_active, 1, MPI_INTEGER, MPI_SUM, flow%comm, ierr)
         if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for active chemistry count')
      end if

      if (global_active <= 0) then
         total_time_local = real(MPI_Wtime(), rk) - t_total_start
         call accumulate_chemistry_screening_stats(nloc, nactive, skipped_temperature_local, &
                                                   skipped_named_species_local, skipped_legacy_mass_fraction_local, &
                                                   total_time_local, reactor_time_local)
         call write_chemistry_diagnostics_row(flow, params, step, dt_local, nloc, nactive, nloc - nactive, &
                                              total_time_local, reactor_time_local, zero, zero, zero, 1.0_rk, &
                                              zero, zero, zero, zero, zero)
         deallocate(active_cell)
         return
      end if

      if (nactive > 0) then
         allocate(T_local(nactive), P_local(nactive), Y_local(params%nspecies, nactive))
         allocate(T_before(nactive), rho_before(nactive), Y_before(params%nspecies, nactive))
         allocate(h_local(nactive), cp_local(nactive), lambda_local(nactive), rho_local(nactive), alpha_local(nactive), qdot_local(nactive))
         alpha_local = 1.0_rk
         P_local = params%background_press

         do i = 1, nactive
            c = active_cell(i)
            T_local(i) = energy%T(c)
            sum_Y = zero
            do k = 1, params%nspecies
               Y_local(k, i) = max(zero, species_Y(k, c))
               sum_Y = sum_Y + Y_local(k, i)
            end do
            if (sum_Y > tiny_safe) then
               Y_local(:, i) = Y_local(:, i) / sum_Y
            else
               Y_local(:, i) = zero
               Y_local(1, i) = 1.0_rk
            end if
            T_before(i) = T_local(i)
            rho_before(i) = max(energy%rho_thermo(c), tiny_safe)
            Y_before(:, i) = Y_local(:, i)
         end do

         call build_c_species_names(params, c_names_flat, n_len)
         call profiler_start('Chemistry_Cantera_ReactorNet')
         t_reactor_start = real(MPI_Wtime(), rk)
         if (params%enable_chemistry_subcycling) then
            dT_max_est = zero
            if (allocated(prev_chem_dT_dt)) then
               do i = 1, nactive
                  c = active_cell(i)
                  dT_max_est = max(dT_max_est, abs(prev_chem_dT_dt(c)) * dt_local)
               end do
            end if

            n_sub = max(1, min(params%max_chemistry_subcycles, &
                               ceiling(dT_max_est / max(params%subcycling_dT_threshold, tiny_safe))))

            dt_sub = dt_local / real(n_sub, rk)
            do j = 1, n_sub
               call cantera_integrate_const_p_chemistry_c(nactive, dt_sub, T_local, P_local, params%nspecies, &
                                                          Y_local, h_local, cp_local, lambda_local, rho_local, &
                                                          params%energy_reference_T, c_names_flat, n_len, &
                                                          params%chemistry_rtol, params%chemistry_atol, &
                                                          params%chemistry_max_steps, merge(1, 0, params%chemistry_energy_enabled))
            end do
         else
            call cantera_integrate_const_p_chemistry_c(nactive, dt_local, T_local, P_local, params%nspecies, &
                                                       Y_local, h_local, cp_local, lambda_local, rho_local, &
                                                       params%energy_reference_T, c_names_flat, n_len, &
                                                       params%chemistry_rtol, params%chemistry_atol, &
                                                       params%chemistry_max_steps, merge(1, 0, params%chemistry_energy_enabled))
         end if
         reactor_time_local = real(MPI_Wtime(), rk) - t_reactor_start
         call profiler_stop('Chemistry_Cantera_ReactorNet')

         raw_max_dT_local = zero
         raw_max_dY_local = zero
         raw_max_rel_drho_local = zero
         do i = 1, nactive
            raw_max_dT_local = max(raw_max_dT_local, abs(T_local(i) - T_before(i)))
            raw_max_dY_local = max(raw_max_dY_local, maxval(abs(Y_local(:, i) - Y_before(:, i))))
            raw_max_rel_drho_local = max(raw_max_rel_drho_local, &
                                         abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe))
         end do

         use_source_limiter = params%chemistry_limit_source_update
         alpha_base = min(1.0_rk, max(zero, params%chemistry_source_relaxation))
         min_alpha_local = 1.0_rk
         qdot_max_local = zero
         qdot_integral_local = zero
         if (use_source_limiter) then
            do i = 1, nactive
               alpha = alpha_base
               dT_abs = abs(T_local(i) - T_before(i))
               if (params%chemistry_max_dT_per_step > zero .and. dT_abs > params%chemistry_max_dT_per_step) then
                  alpha = min(alpha, params%chemistry_max_dT_per_step / max(dT_abs, tiny_safe))
               end if
               dY_abs = maxval(abs(Y_local(:, i) - Y_before(:, i)))
               if (params%chemistry_max_dY_per_step > zero .and. dY_abs > params%chemistry_max_dY_per_step) then
                  alpha = min(alpha, params%chemistry_max_dY_per_step / max(dY_abs, tiny_safe))
               end if
               rel_drho = abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe)
               if (params%chemistry_max_rel_rho_change_per_step > zero .and. &
                   rel_drho > params%chemistry_max_rel_rho_change_per_step) then
                  alpha = min(alpha, params%chemistry_max_rel_rho_change_per_step / max(rel_drho, tiny_safe))
               end if
               alpha = max(zero, min(1.0_rk, alpha))
               T_local(i) = T_before(i) + alpha * (T_local(i) - T_before(i))
               Y_local(:, i) = Y_before(:, i) + alpha * (Y_local(:, i) - Y_before(:, i))
               Y_local(:, i) = max(zero, Y_local(:, i))
               sum_Y = sum(Y_local(:, i))
               if (sum_Y > tiny_safe) then
                  Y_local(:, i) = Y_local(:, i) / sum_Y
               else
                  Y_local(:, i) = zero
                  Y_local(1, i) = 1.0_rk
               end if
               alpha_local(i) = alpha
               min_alpha_local = min(min_alpha_local, alpha)
            end do
            call cantera_update_thermo_c(nactive, T_local, P_local, params%nspecies, Y_local, &
                                         h_local, cp_local, lambda_local, rho_local, &
                                         params%energy_reference_T, c_names_flat, n_len)
         end if

         call cantera_heat_release_rate_c(nactive, T_local, P_local, params%nspecies, Y_local, qdot_local, &
                                          c_names_flat, n_len)

         qdot_max_local = zero
         qdot_integral_local = zero
         do i = 1, nactive
            c = active_cell(i)
            qdot_max_local = max(qdot_max_local, qdot_local(i))
            qdot_integral_local = qdot_integral_local + qdot_local(i) * mesh%cells(c)%volume
         end do

         max_dT_local = zero
         max_dY_local = zero
         max_rel_drho_local = zero
         do i = 1, nactive
            max_dT_local = max(max_dT_local, abs(T_local(i) - T_before(i)))
            max_dY_local = max(max_dY_local, maxval(abs(Y_local(:, i) - Y_before(:, i))))
            max_rel_drho_local = max(max_rel_drho_local, abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe))
         end do

       if (params%enable_chemistry_load_balancing) then
          nfields = 10 + params%nspecies
          allocate(pack_cell(max(1, nactive)))
          allocate(pack_data(nfields, max(1, nactive)))
          
          do i = 1, nactive
             c = active_cell(i)
             pack_cell(i) = c
             pack_data(1, i) = T_local(i)
             pack_data(2, i) = h_local(i)
             pack_data(3, i) = cp_local(i)
             pack_data(4, i) = lambda_local(i)
             pack_data(5, i) = rho_local(i)
             pack_data(6, i) = alpha_local(i)
             pack_data(7, i) = qdot_local(i)
             pack_data(8, i) = T_before(i)
             pack_data(9, i) = rho_before(i)
             pack_data(10, i) = maxval(abs(Y_local(:, i) - Y_before(:, i)))
             pack_data(11:10+params%nspecies, i) = Y_local(1:params%nspecies, i)
          end do
          
          allocate(active_counts(flow%nprocs))
          allocate(active_displs(flow%nprocs))
          allocate(active_matrix_counts(flow%nprocs))
          allocate(active_matrix_displs(flow%nprocs))
          
          call MPI_Allgather(nactive, 1, MPI_INTEGER, &
                             active_counts, 1, MPI_INTEGER, flow%comm, ierr)
          if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgather active counts failed')
          
          active_displs(1) = 0
          active_matrix_counts(1) = active_counts(1) * nfields
          active_matrix_displs(1) = 0
          
          do r = 2, flow%nprocs
             active_displs(r) = active_displs(r - 1) + active_counts(r - 1)
             active_matrix_counts(r) = active_counts(r) * nfields
             active_matrix_displs(r) = active_matrix_displs(r - 1) + active_matrix_counts(r - 1)
          end do
          
          global_active = sum(active_counts)
          
          allocate(all_active_cells(max(1, global_active)))
          allocate(all_data_new(nfields, max(1, global_active)))
          
          ! Gather active cell indices
          call MPI_Allgatherv(pack_cell, nactive, MPI_INTEGER, &
                              all_active_cells, active_counts, active_displs, MPI_INTEGER, flow%comm, ierr)
          if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgatherv active cells failed')
          
          ! Gather packed data
          call MPI_Allgatherv(pack_data, nactive * nfields, MPI_DOUBLE_PRECISION, &
                              all_data_new, active_matrix_counts, active_matrix_displs, MPI_DOUBLE_PRECISION, flow%comm, ierr)
          if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allgatherv active data failed')
          
          ! Unpack gathered data directly into global replicated arrays
          do i = 1, global_active
             c = all_active_cells(i)
             
             energy%T(c) = all_data_new(1, i)
             energy%h(c) = all_data_new(2, i)
             energy%cp(c) = all_data_new(3, i)
             energy%lambda(c) = all_data_new(4, i)
             energy%rho_thermo(c) = all_data_new(5, i)
             
             alpha = all_data_new(6, i)
             qdot_val = all_data_new(7, i)
             T_prev = all_data_new(8, i)
             rho_prev = all_data_new(9, i)
             dY_max_val = all_data_new(10, i)
             
             species_Y(1:params%nspecies, c) = all_data_new(11:10+params%nspecies, i)
             
             if (allocated(energy%chem_active)) energy%chem_active(c) = one
             if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha(c) = alpha
             if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt(c) = (energy%T(c) - T_prev) / max(dt_local, tiny_safe)
             if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt(c) = dY_max_val / max(dt_local, tiny_safe)
             if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt(c) = &
                (abs(energy%rho_thermo(c) - rho_prev) / max(rho_prev, tiny_safe)) / max(dt_local, tiny_safe)
             if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate(c) = qdot_val
          end do
          
          ! Enforce mass fraction normalization
          do i = 1, global_active
             c = all_active_cells(i)
             species_Y(1:params%nspecies, c) = max(zero, species_Y(1:params%nspecies, c))
             sum_Y = sum(species_Y(1:params%nspecies, c))
             if (sum_Y > tiny_safe) then
                species_Y(1:params%nspecies, c) = species_Y(1:params%nspecies, c) / sum_Y
             else
                species_Y(:, c) = zero
                species_Y(1, c) = 1.0_rk
             end if
          end do
          
          deallocate(pack_cell, pack_data)
          deallocate(active_counts, active_displs, active_matrix_counts, active_matrix_displs)
          deallocate(all_active_cells, all_data_new)
       else
          do i = 1, nactive
             c = active_cell(i)
             energy%T(c) = T_local(i)
             energy%h(c) = h_local(i)
             energy%cp(c) = cp_local(i)
             energy%lambda(c) = lambda_local(i)
             energy%rho_thermo(c) = rho_local(i)
             do k = 1, params%nspecies
                species_Y(k, c) = max(zero, Y_local(k, i))
             end do
             sum_Y = sum(species_Y(1:params%nspecies, c))
             if (sum_Y > tiny_safe) then
                species_Y(1:params%nspecies, c) = species_Y(1:params%nspecies, c) / sum_Y
             else
                species_Y(:, c) = zero
                species_Y(1, c) = 1.0_rk
             end if
             if (allocated(energy%chem_active)) energy%chem_active(c) = one
             if (allocated(energy%chem_source_alpha)) energy%chem_source_alpha(c) = alpha_local(i)
             if (allocated(energy%chem_dT_dt)) energy%chem_dT_dt(c) = (T_local(i) - T_before(i)) / max(dt_local, tiny_safe)
             if (allocated(energy%chem_dY_max_dt)) energy%chem_dY_max_dt(c) = &
                maxval(abs(Y_local(:, i) - Y_before(:, i))) / max(dt_local, tiny_safe)
             if (allocated(energy%chem_rel_rho_change_dt)) energy%chem_rel_rho_change_dt(c) = &
                (abs(rho_local(i) - rho_before(i)) / max(rho_before(i), tiny_safe)) / max(dt_local, tiny_safe)
             if (allocated(energy%chem_heat_release_rate)) energy%chem_heat_release_rate(c) = qdot_local(i)
          end do
       end if
    end if
 
    if (.not. allocated(T_local)) then
       max_dT_local = zero
       max_dY_local = zero
       max_rel_drho_local = zero
       raw_max_dT_local = zero
       raw_max_dY_local = zero
       raw_max_rel_drho_local = zero
       min_alpha_local = 1.0_rk
       qdot_max_local = zero
       qdot_integral_local = zero
    end if
    reduce_in = [max_dT_local, max_dY_local, max_rel_drho_local]
    raw_reduce_in = [raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local]
    call MPI_Allreduce(raw_reduce_in, raw_reduce_out, 3, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
    if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for raw chemistry source diagnostics')
    call MPI_Allreduce(reduce_in, reduce_out, 3, MPI_DOUBLE_PRECISION, MPI_MAX, flow%comm, ierr)
    if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for chemistry source diagnostics')
    call MPI_Allreduce(min_alpha_local, alpha_reduce, 1, MPI_DOUBLE_PRECISION, MPI_MIN, flow%comm, ierr)
    if (ierr /= MPI_SUCCESS) call fatal_error('chemistry', 'MPI_Allreduce failed for chemistry source limiter alpha')
    if (present(max_dT_chem)) max_dT_chem = reduce_out(1)
    if (present(max_dY_chem)) max_dY_chem = reduce_out(2)
    if (present(max_rel_drho_chem)) max_rel_drho_chem = reduce_out(3)
    if (present(min_source_alpha)) min_source_alpha = alpha_reduce
    if (present(raw_max_dT_chem)) raw_max_dT_chem = raw_reduce_out(1)
    if (present(raw_max_dY_chem)) raw_max_dY_chem = raw_reduce_out(2)
    if (present(raw_max_rel_drho_chem)) raw_max_rel_drho_chem = raw_reduce_out(3)
 
    if (.not. params%enable_chemistry_load_balancing) then
       call flow_exchange_cell_matrix(flow, species_Y)
       call flow_exchange_cell_scalars(flow, energy%T, energy%h, energy%cp, energy%lambda)
       call flow_exchange_cell_scalar(flow, energy%rho_thermo)
    end if

      total_time_local = real(MPI_Wtime(), rk) - t_total_start
      call accumulate_chemistry_screening_stats(nloc, nactive, skipped_temperature_local, &
                                                skipped_named_species_local, skipped_legacy_mass_fraction_local, &
                                                total_time_local, reactor_time_local, nactive * n_sub)
      call write_chemistry_diagnostics_row(flow, params, step, dt_local, nloc, nactive, nloc - nactive, &
                                           total_time_local, reactor_time_local, &
                                           max_dT_local, max_dY_local, max_rel_drho_local, min_alpha_local, &
                                           raw_max_dT_local, raw_max_dY_local, raw_max_rel_drho_local, &
                                           qdot_max_local, qdot_integral_local)

      if (allocated(T_local)) deallocate(T_local, P_local, Y_local, T_before, Y_before, rho_before, &
                                          h_local, cp_local, lambda_local, rho_local, alpha_local, qdot_local)
      if (allocated(c_names_flat)) deallocate(c_names_flat)
      deallocate(active_cell)
   end subroutine advance_cantera_chemistry