!> app_config から粒子源計画と粒子バッチを構築する実行時変換。
module bem_app_config_particle_runtime
  use bem_kinds, only: dp, i32
  use bem_constants, only: k_boltzmann
  use bem_types, only: mesh_type, particles_soa, sim_config, injection_state, bc_periodic
  use bem_mpi, only: mpi_context, mpi_get_rank_size, mpi_split_count, mpi_bcast_i32_array, mpi_bcast_real_dp_array
  use bem_electrostatic_snapshot, only: electrostatic_snapshot_type
  use bem_collision, only: collision_query_grid_stalled, collision_query_image_limit, &
                           collision_query_index_range, collision_query_invalid_segment, collision_query_ok
  use bem_injection, only: &
    seed_rng, sample_uniform_positions, sample_shifted_maxwell_velocities, compute_macro_particles_for_batch, &
    compute_macro_particles_from_flux, sample_reservoir_face_particles, sample_reservoir_velocity_grid_particles, &
    sample_photo_raycast_particles, &
    compute_inflow_flux_from_drifting_maxwellian, compute_face_area_from_bounds
  use bem_external_boundary_contract, only: &
    external_boundary_contract_type, external_boundary_ok, external_inflow_none, external_inflow_scalar_barrier, &
    resolve_external_boundary_contract
  use bem_app_config_types, only: &
    app_config, particle_species_spec, particles_per_batch_from_config, particle_inflow_reservoir
  use bem_app_config_potential_runtime, only: &
    compute_face_average_potential, warn_face_average_potential_variation, resolve_face_sampling_geometry
  use bem_string_utils, only: lower_ascii
  use bem_config_helpers, only: resolve_inward_normal
  use, intrinsic :: iso_fortran_env, only: error_unit
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  implicit none

  type, public :: particle_source_plan_type
    private
    logical :: ready = .false.
    logical :: mpi_argument_present = .false.
    logical :: use_collective_reservoir_count = .false.
    integer(i32) :: nspecies = 0_i32
    integer(i32) :: mpi_rank = 0_i32
    integer(i32) :: mpi_size = 1_i32
    real(dp), allocatable :: effective_density_m3(:)
    real(dp), allocatable :: effective_particle_flux_m2_s(:)
    logical, allocatable :: number_flux_override_active(:)
    real(dp), allocatable :: effective_temperature_k(:)
    real(dp), allocatable :: effective_drift_velocity(:, :)
    real(dp), allocatable :: effective_weight(:)
    real(dp), allocatable :: photo_emit_current_density(:)
    real(dp), allocatable :: photo_normal_drift_speed(:)
    logical, allocatable :: kinetic_inflow_active(:)
    real(dp), allocatable :: kinetic_reservoir_potential_v(:)
    real(dp), allocatable :: kinetic_access_potential_v(:)
    integer(i32), allocatable :: kinetic_inflow_face(:)
  end type particle_source_plan_type

  interface
    module subroutine init_particle_batch_from_config( &
      cfg, batch_idx, pcls, state, mesh, photo_emission_dq, mpi_rank, mpi_size, mpi, &
      collision_failure_status, collision_failure_species, collision_failure_ray, collision_failure_bounce, snapshot, &
      source_plan, photo_emission_dq_by_species &
      )
      type(app_config), intent(in) :: cfg
      integer(i32), intent(in) :: batch_idx
      type(particles_soa), intent(out) :: pcls
      type(injection_state), intent(inout), optional :: state
      type(mesh_type), intent(in), optional :: mesh
      real(dp), intent(out), optional :: photo_emission_dq(:)
      integer(i32), intent(in), optional :: mpi_rank, mpi_size
      type(mpi_context), intent(in), optional :: mpi
      integer(i32), intent(out), optional :: collision_failure_status, collision_failure_species
      integer(i32), intent(out), optional :: collision_failure_ray, collision_failure_bounce
      type(electrostatic_snapshot_type), intent(inout), optional :: snapshot
      type(particle_source_plan_type), intent(in), optional, target :: source_plan
      real(dp), intent(out), optional :: photo_emission_dq_by_species(:, :)
    end subroutine init_particle_batch_from_config

    module subroutine sample_species_state( &
      sim, spec, n, x, v, barrier_normal_energy, vmin_normal, apply_barrier_energy_shift, &
      temperature_k_override, drift_velocity_override &
      )
      type(sim_config), intent(in) :: sim
      type(particle_species_spec), intent(in) :: spec
      integer(i32), intent(in) :: n
      real(dp), intent(out) :: x(:, :)
      real(dp), intent(out) :: v(:, :)
      real(dp), intent(in), optional :: barrier_normal_energy
      real(dp), intent(in), optional :: vmin_normal
      logical, intent(in), optional :: apply_barrier_energy_shift
      real(dp), intent(in), optional :: temperature_k_override
      real(dp), intent(in), optional :: drift_velocity_override(3)
    end subroutine sample_species_state

    module subroutine sample_photo_species_state( &
      sim, spec, mesh, n_rays, x, v, w, n_emit, emit_elem_idx, global_rays_per_batch, &
      emit_current_density_override, normal_drift_speed_override, vmin_normal, &
      collision_failure_status, collision_failure_ray, collision_failure_bounce &
      )
      type(sim_config), intent(in) :: sim
      type(particle_species_spec), intent(in) :: spec
      type(mesh_type), intent(in) :: mesh
      integer(i32), intent(in) :: n_rays
      real(dp), intent(out) :: x(:, :)
      real(dp), intent(out) :: v(:, :)
      real(dp), intent(out) :: w(:)
      integer(i32), intent(out) :: n_emit
      integer(i32), intent(out), optional :: emit_elem_idx(:)
      integer(i32), intent(in), optional :: global_rays_per_batch
      real(dp), intent(in), optional :: emit_current_density_override
      real(dp), intent(in), optional :: normal_drift_speed_override
      real(dp), intent(in), optional :: vmin_normal
      integer(i32), intent(out), optional :: collision_failure_status, collision_failure_ray, collision_failure_bounce
    end subroutine sample_photo_species_state

    module subroutine reservoir_face_velocity_correction( &
      cfg, spec, vmin_normal, barrier_normal, mesh, snapshot, warn_face_variation, boundary_contract &
      )
      type(app_config), intent(in) :: cfg
      type(particle_species_spec), intent(in) :: spec
      real(dp), intent(out) :: vmin_normal
      real(dp), intent(out) :: barrier_normal
      type(mesh_type), intent(in), optional :: mesh
      type(electrostatic_snapshot_type), intent(inout), optional :: snapshot
      logical, intent(in), optional :: warn_face_variation
      type(external_boundary_contract_type), intent(in), optional :: boundary_contract
    end subroutine reservoir_face_velocity_correction

    module subroutine external_kinetic_face_velocity_correction( &
      cfg, spec, reservoir_potential_v, access_potential_v, vmin_normal, barrier_normal, &
      mesh, snapshot, warn_face_variation &
      )
      type(app_config), intent(in) :: cfg
      type(particle_species_spec), intent(in) :: spec
      real(dp), intent(in) :: reservoir_potential_v, access_potential_v
      real(dp), intent(out) :: vmin_normal, barrier_normal
      type(mesh_type), intent(in), optional :: mesh
      type(electrostatic_snapshot_type), intent(inout), optional :: snapshot
      logical, intent(in), optional :: warn_face_variation
    end subroutine external_kinetic_face_velocity_correction
  end interface

contains

  !> バッチ生成前に乱数シードだけを初期化する。
  !! @param[in] cfg 乱数シード値 `sim.rng_seed` を含むアプリ設定。
  subroutine seed_particles_from_config(cfg, mpi_rank, mpi_size, mpi)
    type(app_config), intent(in) :: cfg
    integer(i32), intent(in), optional :: mpi_rank, mpi_size
    type(mpi_context), intent(in), optional :: mpi
    integer(i32) :: local_rank, n_ranks, seed_value
    integer(kind=8) :: seed_tmp

    call resolve_parallel_rank_size(local_rank, n_ranks, mpi_rank, mpi_size, mpi, 'seed_particles_from_config')

    seed_tmp = int(cfg%sim%rng_seed, kind=8) + 104729_8*int(local_rank, kind=8)
    seed_value = int(modulo(seed_tmp, int(huge(0_i32), kind=8)), kind=i32)
    call seed_rng([seed_value])
  end subroutine seed_particles_from_config

  !> 設定とMPI配置だけに依存する粒子 source の導出値を構築する。
  !! 乱数、残差、mesh/snapshot依存の障壁は扱わず、run中に不変な値だけを保持する。
  subroutine build_particle_source_plan( &
    cfg, plan, mpi_rank, mpi_size, mpi, kinetic_inflow_active, kinetic_reservoir_potential_v, &
    kinetic_access_potential_v, kinetic_inflow_face, number_flux_override_active, number_flux_override_m2_s &
    )
    type(app_config), intent(in) :: cfg
    type(particle_source_plan_type), intent(out) :: plan
    integer(i32), intent(in), optional :: mpi_rank, mpi_size
    type(mpi_context), intent(in), optional :: mpi
    logical, intent(in), optional :: kinetic_inflow_active(:)
    real(dp), intent(in), optional :: kinetic_reservoir_potential_v(:), kinetic_access_potential_v(:)
    integer(i32), intent(in), optional :: kinetic_inflow_face(:)
    logical, intent(in), optional :: number_flux_override_active(:)
    real(dp), intent(in), optional :: number_flux_override_m2_s(:)

    integer(i32) :: s, local_rank, n_ranks
    logical :: has_enabled_reservoir

    call resolve_parallel_rank_size(local_rank, n_ranks, mpi_rank, mpi_size, mpi, 'build_particle_source_plan')
    plan%nspecies = cfg%n_particle_species
    plan%mpi_rank = local_rank
    plan%mpi_size = n_ranks
    plan%mpi_argument_present = present(mpi)
    has_enabled_reservoir = .false.
    do s = 1, cfg%n_particle_species
      if (.not. cfg%particle_species(s)%enabled) cycle
      has_enabled_reservoir = has_enabled_reservoir .or. &
                              trim(lower_ascii(cfg%particle_species(s)%source_mode)) == 'reservoir_face' .or. &
                              trim(lower_ascii(cfg%particle_species(s)%source_mode)) == 'plane_source' .or. &
                              has_boundary_inflow(cfg%particle_species(s))
    end do
    plan%use_collective_reservoir_count = present(mpi) .and. has_enabled_reservoir

    allocate (plan%effective_density_m3(cfg%n_particle_species))
    allocate (plan%effective_particle_flux_m2_s(cfg%n_particle_species))
    allocate (plan%number_flux_override_active(cfg%n_particle_species))
    allocate (plan%effective_temperature_k(cfg%n_particle_species))
    allocate (plan%effective_drift_velocity(3, cfg%n_particle_species))
    allocate (plan%effective_weight(cfg%n_particle_species))
    allocate (plan%photo_emit_current_density(cfg%n_particle_species))
    allocate (plan%photo_normal_drift_speed(cfg%n_particle_species))
    allocate (plan%kinetic_inflow_active(cfg%n_particle_species))
    allocate (plan%kinetic_reservoir_potential_v(cfg%n_particle_species))
    allocate (plan%kinetic_access_potential_v(cfg%n_particle_species))
    allocate (plan%kinetic_inflow_face(cfg%n_particle_species))
    plan%effective_density_m3 = 0.0_dp
    plan%effective_particle_flux_m2_s = 0.0_dp
    plan%number_flux_override_active = .false.
    plan%effective_temperature_k = 0.0_dp
    plan%effective_drift_velocity = 0.0_dp
    plan%effective_weight = 0.0_dp
    plan%photo_emit_current_density = 0.0_dp
    plan%photo_normal_drift_speed = 0.0_dp
    plan%kinetic_inflow_active = .false.
    plan%kinetic_reservoir_potential_v = 0.0_dp
    plan%kinetic_access_potential_v = 0.0_dp
    plan%kinetic_inflow_face = 0_i32

    if (present(kinetic_inflow_active) .or. present(kinetic_reservoir_potential_v) .or. &
        present(kinetic_access_potential_v) .or. present(kinetic_inflow_face)) then
      if (.not. present(kinetic_inflow_active) .or. .not. present(kinetic_reservoir_potential_v) .or. &
          .not. present(kinetic_access_potential_v) .or. .not. present(kinetic_inflow_face)) then
        error stop 'particle source kinetic map requires active, reservoir, access, and face arrays together.'
      end if
      if (size(kinetic_inflow_active) /= cfg%n_particle_species .or. &
          size(kinetic_reservoir_potential_v) /= cfg%n_particle_species .or. &
          size(kinetic_access_potential_v) /= cfg%n_particle_species .or. &
          size(kinetic_inflow_face) /= cfg%n_particle_species) then
        error stop 'particle source kinetic map species count mismatch.'
      end if
      if (.not. all(ieee_is_finite(kinetic_reservoir_potential_v)) .or. &
          .not. all(ieee_is_finite(kinetic_access_potential_v))) then
        error stop 'particle source kinetic map potentials must be finite.'
      end if
      plan%kinetic_inflow_active = kinetic_inflow_active
      plan%kinetic_reservoir_potential_v = kinetic_reservoir_potential_v
      plan%kinetic_access_potential_v = kinetic_access_potential_v
      plan%kinetic_inflow_face = kinetic_inflow_face
      if (any(plan%kinetic_inflow_active .and. &
              (plan%kinetic_inflow_face < 1_i32 .or. plan%kinetic_inflow_face > 6_i32))) then
        error stop 'active particle source kinetic map face must be in [1, 6].'
      end if
    end if

    if (present(number_flux_override_active) .or. present(number_flux_override_m2_s)) then
      if (.not. present(number_flux_override_active) .or. .not. present(number_flux_override_m2_s)) then
        error stop 'particle source number-flux override requires active and flux arrays together.'
      end if
      if (size(number_flux_override_active) /= cfg%n_particle_species .or. &
          size(number_flux_override_m2_s) /= cfg%n_particle_species) then
        error stop 'particle source number-flux override species count mismatch.'
      end if
      if (any(number_flux_override_active .and. &
              (.not. ieee_is_finite(number_flux_override_m2_s) .or. number_flux_override_m2_s < 0.0_dp))) then
        error stop 'active particle source number-flux overrides must be finite and nonnegative.'
      end if
      plan%number_flux_override_active = number_flux_override_active
    end if

    do s = 1, cfg%n_particle_species
      if (.not. cfg%particle_species(s)%enabled) cycle
      select case (trim(lower_ascii(cfg%particle_species(s)%source_mode)))
      case ('volume_seed')
        plan%effective_weight(s) = cfg%particle_species(s)%w_particle
        plan%effective_temperature_k(s) = species_temperature_k(cfg%particle_species(s))
        plan%effective_drift_velocity(:, s) = cfg%particle_species(s)%drift_velocity
        if (has_boundary_inflow(cfg%particle_species(s))) then
          if (trim(lower_ascii(cfg%particle_species(s)%velocity_distribution)) == 'grid') then
            plan%effective_particle_flux_m2_s(s) = cfg%particle_species(s)%particle_flux_m2_s
          else
            plan%effective_density_m3(s) = species_number_density_m3(cfg%particle_species(s))
          end if
        end if
      case ('reservoir_face', 'plane_source')
        if (trim(lower_ascii(cfg%particle_species(s)%velocity_distribution)) == 'grid') then
          plan%effective_particle_flux_m2_s(s) = cfg%particle_species(s)%particle_flux_m2_s
        else
          plan%effective_density_m3(s) = species_number_density_m3(cfg%particle_species(s))
        end if
        plan%effective_weight(s) = cfg%particle_species(s)%w_particle
        plan%effective_temperature_k(s) = species_temperature_k(cfg%particle_species(s))
        plan%effective_drift_velocity(:, s) = cfg%particle_species(s)%drift_velocity
      case ('photo_raycast')
        plan%photo_emit_current_density(s) = cfg%particle_species(s)%emit_current_density_a_m2
        plan%photo_normal_drift_speed(s) = cfg%particle_species(s)%normal_drift_speed
      end select
      if (plan%number_flux_override_active(s)) then
        plan%effective_particle_flux_m2_s(s) = number_flux_override_m2_s(s)
      end if
    end do

    plan%ready = .true.
  end subroutine build_particle_source_plan

  !> reservoir/plane sourceの初期位置を、設定されたbox境界条件に従って有効領域へ正規化する。
  pure subroutine normalize_reservoir_positions(sim, x)
    type(sim_config), intent(in) :: sim
    real(dp), intent(inout) :: x(:, :)
    integer(i32) :: axis
    real(dp) :: span

    if (.not. sim%use_box) return
    do axis = 1_i32, 3_i32
      if (.not. ieee_is_finite(sim%box_min(axis)) .or. .not. ieee_is_finite(sim%box_max(axis))) cycle
      span = sim%box_max(axis) - sim%box_min(axis)
      if (.not. ieee_is_finite(span) .or. span <= 0.0_dp) cycle
      if (sim%bc_low(axis) == bc_periodic .and. sim%bc_high(axis) == bc_periodic) then
        x(axis, :) = sim%box_min(axis) + modulo(x(axis, :) - sim%box_min(axis), span)
      else
        x(axis, :) = min(max(x(axis, :), sim%box_min(axis)), sim%box_max(axis))
      end if
    end do
  end subroutine normalize_reservoir_positions

  !> speciesに外部 reservoir 流入を指定したbox面があるかを返す。
  pure logical function has_boundary_inflow(spec) result(has_inflow)
    type(particle_species_spec), intent(in) :: spec

    has_inflow = any(spec%boundary_inflow_low == particle_inflow_reservoir) .or. &
                 any(spec%boundary_inflow_high == particle_inflow_reservoir)
  end function has_boundary_inflow

  !> face bit順の面に reservoir 流入が有効かを返す。
  pure logical function boundary_inflow_face_enabled(spec, face) result(enabled)
    type(particle_species_spec), intent(in) :: spec
    integer(i32), intent(in) :: face

    if (mod(face, 2_i32) == 1_i32) then
      enabled = spec%boundary_inflow_low((face + 1_i32)/2_i32) == particle_inflow_reservoir
    else
      enabled = spec%boundary_inflow_high(face/2_i32) == particle_inflow_reservoir
    end if
  end function boundary_inflow_face_enabled

  !> box全面を開口とする一時的な legacy reservoir spec を構築する。
  pure subroutine make_boundary_inflow_spec(sim, source_spec, face, inflow_spec)
    type(sim_config), intent(in) :: sim
    type(particle_species_spec), intent(in) :: source_spec
    integer(i32), intent(in) :: face
    type(particle_species_spec), intent(out) :: inflow_spec
    integer(i32) :: axis
    logical :: high_side

    inflow_spec = source_spec
    inflow_spec%source_mode = 'reservoir_face'
    call boundary_face_name(face, inflow_spec%inject_face)
    inflow_spec%pos_low = sim%box_min
    inflow_spec%pos_high = sim%box_max
    axis = (face + 1_i32)/2_i32
    high_side = mod(face, 2_i32) == 0_i32
    if (high_side) then
      inflow_spec%pos_low(axis) = sim%box_max(axis)
      inflow_spec%pos_high(axis) = sim%box_max(axis)
    else
      inflow_spec%pos_low(axis) = sim%box_min(axis)
      inflow_spec%pos_high(axis) = sim%box_min(axis)
    end if
    inflow_spec%boundary_inflow_low = 0_i32
    inflow_spec%boundary_inflow_high = 0_i32
  end subroutine make_boundary_inflow_spec

  !> plane_sourceの内部面をlegacy face samplerの仮想box境界へ写像する。
  pure subroutine configure_plane_source_box(spec, box_min, box_max)
    type(particle_species_spec), intent(in) :: spec
    real(dp), intent(inout) :: box_min(3), box_max(3)
    integer(i32) :: axis

    select case (trim(lower_ascii(spec%inject_face)))
    case ('x_low')
      axis = 1_i32
      box_min(axis) = spec%pos_low(axis)
    case ('x_high')
      axis = 1_i32
      box_max(axis) = spec%pos_low(axis)
    case ('y_low')
      axis = 2_i32
      box_min(axis) = spec%pos_low(axis)
    case ('y_high')
      axis = 2_i32
      box_max(axis) = spec%pos_low(axis)
    case ('z_low')
      axis = 3_i32
      box_min(axis) = spec%pos_low(axis)
    case ('z_high')
      axis = 3_i32
      box_max(axis) = spec%pos_low(axis)
    case default
      error stop 'plane_source has invalid derived face.'
    end select
  end subroutine configure_plane_source_box

  !> face bit順の面名を返す。
  pure subroutine boundary_face_name(face, name)
    integer(i32), intent(in) :: face
    character(len=*), intent(out) :: name

    select case (face)
    case (1_i32)
      name = 'x_low'
    case (2_i32)
      name = 'x_high'
    case (3_i32)
      name = 'y_low'
    case (4_i32)
      name = 'y_high'
    case (5_i32)
      name = 'z_low'
    case (6_i32)
      name = 'z_high'
    case default
      error stop 'invalid boundary inflow face.'
    end select
  end subroutine boundary_face_name

  !> 面名をboundary bit順のindexへ変換する。
  pure integer(i32) function injection_face_index(face_name) result(face)
    character(len=*), intent(in) :: face_name

    select case (trim(lower_ascii(face_name)))
    case ('x_low')
      face = 1_i32
    case ('x_high')
      face = 2_i32
    case ('y_low')
      face = 3_i32
    case ('y_high')
      face = 4_i32
    case ('z_low')
      face = 5_i32
    case ('z_high')
      face = 6_i32
    case default
      face = 0_i32
    end select
  end function injection_face_index

  !> reservoir_face 用に、物理流量と残差から今バッチのマクロ粒子数を決める。
  !! @param[in] sim ボックス境界・バッチ時間などのシミュレーション設定。
  !! @param[in] spec reservoir_face 粒子種設定。
  !! @param[inout] residual 前バッチから繰り越した端数。
  !! @param[out] count 今バッチで生成するマクロ粒子数。
  !! @param[in] vmin_normal 法線速度の下限 [m/s]（省略時は 0）。
  !! @param[in] number_density_override 数密度の上書き値 [1/m^3]。
  !! @param[in] particle_flux_override 粒子数 flux の上書き値 [1/m^2/s]。
  !! @param[in] w_particle_override マクロ粒子重みの上書き値。
  !! @param[in] temperature_k_override 温度の上書き値 [K]。
  !! @param[in] drift_velocity_override ドリフト速度の上書き値 [m/s]。
  subroutine compute_macro_particles_for_species( &
    sim, spec, residual, count, vmin_normal, number_density_override, w_particle_override, &
    temperature_k_override, drift_velocity_override, particle_flux_override, use_particle_flux_override &
    )
    type(sim_config), intent(in) :: sim
    type(particle_species_spec), intent(in) :: spec
    real(dp), intent(inout) :: residual
    integer(i32), intent(out) :: count
    real(dp), intent(in), optional :: vmin_normal
    real(dp), intent(in), optional :: number_density_override
    real(dp), intent(in), optional :: w_particle_override
    real(dp), intent(in), optional :: temperature_k_override
    real(dp), intent(in), optional :: drift_velocity_override(3)
    real(dp), intent(in), optional :: particle_flux_override
    logical, intent(in), optional :: use_particle_flux_override

    real(dp) :: number_density_m3, effective_batch_duration, particle_flux_m2_s, w_particle, temperature_k_local
    real(dp) :: drift_velocity_local(3)
    logical :: direct_flux

    number_density_m3 = species_number_density_m3(spec)
    if (present(number_density_override)) number_density_m3 = number_density_override
    w_particle = spec%w_particle
    if (present(w_particle_override)) w_particle = w_particle_override
    temperature_k_local = species_temperature_k(spec)
    if (present(temperature_k_override)) temperature_k_local = temperature_k_override
    drift_velocity_local = spec%drift_velocity
    if (present(drift_velocity_override)) drift_velocity_local = drift_velocity_override
    effective_batch_duration = sim%batch_duration
    direct_flux = trim(lower_ascii(spec%velocity_distribution)) == 'grid'
    if (present(use_particle_flux_override)) direct_flux = direct_flux .or. use_particle_flux_override
    if (direct_flux) then
      particle_flux_m2_s = spec%particle_flux_m2_s
      if (present(particle_flux_override)) particle_flux_m2_s = particle_flux_override
      call compute_macro_particles_from_flux( &
        particle_flux_m2_s, spec%inject_face, spec%pos_low, spec%pos_high, effective_batch_duration, w_particle, residual, count &
        )
      return
    end if
    if (present(vmin_normal)) then
      call compute_macro_particles_for_batch( &
        number_density_m3, temperature_k_local, spec%m_particle, drift_velocity_local, sim%box_min, sim%box_max, &
        spec%inject_face, spec%pos_low, spec%pos_high, effective_batch_duration, w_particle, residual, count, &
        vmin_normal=vmin_normal &
        )
    else
      call compute_macro_particles_for_batch( &
        number_density_m3, temperature_k_local, spec%m_particle, drift_velocity_local, sim%box_min, sim%box_max, &
        spec%inject_face, spec%pos_low, spec%pos_high, effective_batch_duration, w_particle, residual, count &
        )
    end if
  end subroutine compute_macro_particles_for_species

  !> 粒子種設定から実効密度[m^-3]を返す。
  !! @param[in] spec 粒子種設定。
  !! @return number_density_m3 実効粒子数密度 [1/m^3]。
  pure real(dp) function species_number_density_m3(spec) result(number_density_m3)
    type(particle_species_spec), intent(in) :: spec

    number_density_m3 = spec%number_density_m3
    if (spec%has_number_density_cm3) number_density_m3 = spec%number_density_cm3*1.0d6
  end function species_number_density_m3

  !> 併存対応のため `mpi_context` と rank/size の両方を受け、最終的なrank/sizeを解決する。
  subroutine resolve_parallel_rank_size(local_rank, n_ranks, mpi_rank, mpi_size, mpi, caller_name)
    integer(i32), intent(out) :: local_rank, n_ranks
    integer(i32), intent(in), optional :: mpi_rank, mpi_size
    type(mpi_context), intent(in), optional :: mpi
    character(len=*), intent(in) :: caller_name

    call mpi_get_rank_size(local_rank, n_ranks, mpi)
    if (present(mpi_rank)) local_rank = mpi_rank
    if (present(mpi_size)) n_ranks = mpi_size
    if (n_ranks <= 0_i32) error stop 'mpi_size must be > 0 in '//trim(caller_name)//'.'
    if (local_rank < 0_i32 .or. local_rank >= n_ranks) then
      error stop 'mpi_rank is out of range in '//trim(caller_name)//'.'
    end if
  end subroutine resolve_parallel_rank_size

  !> 粒子種設定から実効温度[K]を返す。
  !! @param[in] spec 粒子種設定。
  !! @return temperature_k 実効温度 [K]。
  pure real(dp) function species_temperature_k(spec) result(temperature_k)
    type(particle_species_spec), intent(in) :: spec

    temperature_k = spec%temperature_k
    if (spec%has_temperature_ev) temperature_k = spec%temperature_ev*1.160451812d4
  end function species_temperature_k

end module bem_app_config_particle_runtime