advance_species_transport Subroutine

   subroutine advance_species_transport(mesh, flow, bc, params, fields, species, transport)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(bc_set_t), intent(in) :: bc
      type(case_params_t), intent(in) :: params
      type(flow_fields_t), intent(in) :: fields
      type(species_fields_t), intent(inout) :: species
      type(transport_properties_t), intent(in) :: transport

      real(rk), allocatable :: dY(:), diff_flux(:), adv_flux(:), Y_face_lin(:)
      real(rk) :: flux, face_area, dist, sum_diff_flux
      real(rk) :: Y_cell, Y_other, Y_face
      real(rk) :: diff_face
      integer :: c, f, fid, neighbor
      integer :: k
      real(rk) :: sum_Y
      logical :: is_dirichlet

      if (species%nspecies <= 0) return

      allocate(dY(species%nspecies))
      allocate(diff_flux(species%nspecies))
      allocate(adv_flux(species%nspecies))
      allocate(Y_face_lin(species%nspecies))

      species%Y_old = species%Y

      ! Iterate through MPI-owned cells
      do c = flow%first_cell, flow%last_cell
         dY = zero

         do f = 1, mesh%ncell_faces(c)
            fid = mesh%cell_faces(f,c)
            if (mesh%faces(fid)%owner == c) then
               flux = fields%face_flux(fid)
            else
               flux = -fields%face_flux(fid)
            end if
            neighbor = face_effective_neighbor(mesh, bc, fid, c)

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

            sum_diff_flux = zero
            do k = 1, species%nspecies
               Y_cell = species%Y_old(k, c)
               
               if (neighbor == 0) then
                  call boundary_species(mesh, bc, fid, k, Y_cell, Y_other, is_dirichlet)
               else
                  Y_other = species%Y_old(k, neighbor)
                  is_dirichlet = .true.
               end if

               ! 1. Advective flux using Upwind discretization.
               ! flux is oriented outward from the current cell.
               if (flux > zero) then
                  Y_face = Y_cell
               else
                  Y_face = Y_other
               end if
               adv_flux(k) = -flux * Y_face

               ! 2. Diffusive flux using central difference
               diff_flux(k) = zero
               if (neighbor /= 0 .or. is_dirichlet) then
                  if (neighbor == 0) then
                     diff_face = transport%diffusivity(k, c)
                  else
                     diff_face = 0.5_rk * (transport%diffusivity(k, c) + transport%diffusivity(k, neighbor))
                  end if
                  diff_flux(k) = diff_face * (Y_other - Y_cell) / dist * face_area
               end if
               
               sum_diff_flux = sum_diff_flux + diff_flux(k)
               Y_face_lin(k) = 0.5_rk * (Y_cell + Y_other)
            end do

            ! 3. Apply Correction Velocity to ensure mass conservation: sum(j_k) = 0
            do k = 1, species%nspecies
               dY(k) = dY(k) + adv_flux(k) + (diff_flux(k) - Y_face_lin(k) * sum_diff_flux)
            end do
         end do

         ! Explicit timestep update
         do k = 1, species%nspecies
            species%Y(k,c) = species%Y_old(k,c) + params%dt * dY(k) / mesh%cells(c)%volume
            
            ! Ensure boundedness: Y_k must remain in [0, 1]
            if (species%Y(k,c) < zero) species%Y(k,c) = zero
            if (species%Y(k,c) > one)  species%Y(k,c) = one
         end do

         ! Local renormalization: sum(Y_k) = 1
         sum_Y = zero
         do k = 1, species%nspecies
            sum_Y = sum_Y + species%Y(k,c)
         end do

         if (sum_Y > zero) then
            do k = 1, species%nspecies
               species%Y(k,c) = species%Y(k,c) / sum_Y
            end do
         end if
      end do

      deallocate(dY)
      deallocate(diff_flux)
      deallocate(adv_flux)
      deallocate(Y_face_lin)

      ! Synchronize updated species ghosts for the next transport/property step.
      call flow_exchange_cell_matrix(flow, species%Y)

   end subroutine advance_species_transport