!> Parsing and validation of simulation parameters from `case.nml`.
!!
!! This module handles the ingestion of user-defined parameters using Fortran 
!! namelists. It is structured to read parameters in logical groups:
!! - **`mesh_input`**: Path to the native hexahedral grid files.
!! - **`time_input`**: Step counts, timestep size, and CFL control.
!! - **`fluid_input`**: Constant/Cantera properties (density, viscosity).
!! - **`solver_input`**: Linear solver tolerances and numerical schemes.
!! - **`boundary_input`**: Mapping of BC types and values to mesh patches.
!! - **`species_input`**: Multi-species transport and chemical mechanism selection.
!! - **`output_input`**: Control over VTK and diagnostic file generation.
!! - **`profiling_input`**: Control over profiling enable/disable and nested reporting.
module mod_input
   use mod_precision, only : rk, zero, one, name_len, path_len, fatal_error, lowercase
   implicit none

   private

   integer, parameter, public :: max_patches = 64  !< Maximum number of boundary patches supported.
   integer, parameter, public :: max_species = 256 !< Maximum number of chemical species supported.

   !> Global container for all parsed simulation parameters.
   !!
   !! This structure acts as the source of truth for the entire solver. 
   !! It is populated once at startup and accessed as a read-only object 
   !! by the various solver modules.
   type, public :: case_params_t
      !> @name Mesh Configuration
      character(len=path_len) :: mesh_dir = "mesh_native" !< Directory containing `points.dat`, `cells.dat`, etc.

      !> @name Time Stepping
      integer :: nsteps = 0              !< Total number of timesteps to execute.
      real(rk) :: dt = zero              !< Fixed timestep size [s].
      real(rk) :: dt_old = zero          !< Timestep size used in the previous step [s].
      integer :: output_interval = 1     !< Frequency of VTK/PVD output (in steps).
      integer :: terminal_interval = 1   !< Frequency of lightweight terminal status updates (in steps).
      character(len=name_len) :: terminal_detail = "health" !< Terminal verbosity: brief, health, or full.
      logical :: use_dynamic_dt = .false.!< If true, `dt` scales to maintain target `max_cfl`.
      real(rk) :: max_cfl = 0.5_rk       !< Target maximum CFL number for stability.
      real(rk) :: max_dt = huge(one)     !< Absolute timestep ceiling for dynamic dt [s]; huge disables.
      real(rk) :: min_dt = 1.0e-12_rk    !< Absolute timestep floor for automatic dt cuts [s].
      real(rk) :: dt_growth_limit = 1.02_rk !< Per-step dynamic-dt growth factor ceiling.
      logical :: enable_step_rejection = .false. !< If true, reject and retry a timestep when health checks fail.
      integer :: max_step_retries = 6 !< Maximum rejected attempts for a single timestep.
      real(rk) :: step_reject_cut_factor = 0.5_rk !< Multiplicative dt cut applied after a rejected trial step.
      real(rk) :: step_reject_min_dt = 1.0e-12_rk !< Absolute lower bound for rejection retries [s].
      real(rk) :: max_KE_growth_per_step = 10.0_rk !< Reject if KE grows by more than this factor; <=0 disables.
      real(rk) :: ke_reject_floor = 1.0e-12_rk !< Do not apply relative KE-growth rejection when KE before the step is below this floor.
      real(rk) :: max_cfl_overshoot_factor = 1.25_rk !< Reject if post-step CFL exceeds factor*max_cfl; <=0 disables.
      logical :: reject_on_pressure_max_iter = .true. !< Reject if the pressure solve reaches pressure_max_iter.
      logical :: stop_on_min_dt_failure = .true. !< Stop if a rejected step can no longer reduce dt.

      !> @name Fluid Properties
      real(rk) :: rho = one              !< Constant active flow/projection density [kg/m^3] when `enable_variable_density=.false.`.
      logical :: enable_variable_density = .false. !< Enable experimental variable-density low-Mach coupling.
      logical :: enable_density_corrector_projection = .false. !< Post-chemistry pressure-only projection to re-enforce div(u)=S after density change.
      logical :: upwind_dilatation_density = .false. !< Use upwind density in the advective divergence source calculation.
      character(len=name_len) :: density_eos = "constant" !< Density source selector; `cantera` uses the selected YAML phase density.
      real(rk) :: nu = 1.0e-2_rk         !< Constant kinematic viscosity [m^2/s].
      logical :: enable_variable_nu = .false.    !< Allow Cantera-updated flow viscosity/nu; default keeps validation Re fixed.
      integer :: transport_update_interval = 1 !< Cantera transport-property update interval for mu/D_k only [steps].

      !> @name Linear Solver & Numerics
      integer :: pressure_max_iter = 300 !< Maximum Conjugate Gradient iterations for the Poisson solver.
      real(rk) :: pressure_tol = 1.0e-10_rk !< Relative residual tolerance for pressure convergence.
      real(rk) :: pressure_abs_tol = 1.0e-8_rk !< Absolute RMS residual tolerance; PCG exits when either criterion is met.
      real(rk) :: body_force(3) = zero   !< Constant volumetric acceleration vector \((a_x, a_y, a_z)\) [m/s^2].
      character(len=name_len) :: convection_scheme = "upwind" !< Legacy momentum advection selector: "upwind" or "central".
      character(len=name_len) :: momentum_convection_scheme = "" !< Optional momentum-specific scheme; blank falls back to convection_scheme.
      character(len=name_len) :: species_convection_scheme = "upwind" !< Species advection: upwind, central, bounded_central, bounded_linear, muscl, limited_linear.
      character(len=name_len) :: energy_convection_scheme = "upwind" !< Enthalpy advection: upwind, central, bounded_central, bounded_linear, muscl, limited_linear.
      character(len=name_len) :: scalar_limiter = "barth_jespersen" !< Local-bounds limiter for high-order scalar reconstruction.
      character(len=name_len) :: species_time_scheme = "euler" !< Species transport time scheme: euler or ab2.
      character(len=name_len) :: energy_time_scheme = "euler" !< Enthalpy transport time scheme: euler or ab2.

      !> @name Boundary Condition Mapping
      integer :: n_patches = 0           !< Total number of patches defined in the namelist.
      character(len=name_len) :: patch_name(max_patches) = "" !< Names identifying mesh patches.
      character(len=name_len) :: patch_type(max_patches) = "" !< Legacy BC type string.

      !> Field-specific BC Overrides
      character(len=name_len) :: patch_velocity_type(max_patches) = "" !< Override BC type for velocity.
      character(len=name_len) :: patch_pressure_type(max_patches) = "" !< Override BC type for pressure.
      character(len=name_len) :: patch_temperature_type(max_patches) = "" !< Override BC type for temperature/enthalpy.

      real(rk) :: patch_u(max_patches) = zero    !< Specified x-velocity on patch [m/s].
      real(rk) :: patch_v(max_patches) = zero    !< Specified y-velocity on patch [m/s].
      real(rk) :: patch_w(max_patches) = zero    !< Specified z-velocity on patch [m/s].
      !< Mass flux input [kg/m^2/s]. Signed for mass_flux;
      !! positive magnitude for inlet/outlet_mass_flux.
      real(rk) :: patch_mass_flux(max_patches) = zero
      real(rk) :: patch_p(max_patches) = zero    !< Specified projection-pressure boundary value [Pa], not Cantera \(p_0\).
      real(rk) :: patch_dpdn(max_patches) = zero !< Specified pressure gradient on patch [Pa/m].
      real(rk) :: patch_T(max_patches) = 300.0_rk !< Specified temperature on patch [K].

      !> Species Boundary Conditions
      character(len=name_len) :: patch_species_type(max_patches) = "" !< Override BC type for mass fractions.
      real(rk) :: patch_Y(max_species, max_patches) = zero           !< Legacy index-based \(Y_k\) mass fractions on patch.
      character(len=1024) :: patch_composition(max_patches) = ""      !< Preferred named patch mixture, e.g. "CH4:1,O2:0.233,N2:0.767".

      !> @name Data Output
      character(len=path_len) :: output_dir = "output" !< Directory for result storage.
      logical :: write_vtu = .true.                    !< If true, generates Unstructured VTK files.
      logical :: write_diagnostics = .true.            !< If true, writes global residuals to `diagnostics.csv`.
      character(len=name_len) :: vtu_format = "ascii"  !< Output format: "ascii" or "binary".

      !> @name Restart Controls
      logical :: restart_from_file = .false.           !< If true, initialize fields from restart_file after allocation.
      character(len=path_len) :: restart_file = ""     !< Rank-independent restart file to read.
      logical :: write_restart = .false.               !< If true, periodically write global restart files.
      integer :: restart_interval = 0                  !< Restart write cadence [steps]; required when write_restart=.true.
      character(len=path_len) :: restart_output_dir = "restart" !< Restart directory, relative to output_dir unless absolute.
      integer :: restart_keep = 0                      !< Number of newest restart files to keep; 0 keeps all.
      logical :: restart_use_file_dt = .false.         !< If true, restore dt from restart; otherwise namelist dt wins.
      logical :: restart_nsteps_are_additional = .false. !< If true, add namelist nsteps to restart step.

      !> @name Multi-Species Transport
      logical :: enable_species = .false.                         !< Enable advection-diffusion of mass fractions.
      logical :: enable_reactions = .false.                       !< Reserved for chemical source terms; supported paths require non-reacting transport.
      integer :: nspecies = 0                                     !< Total number of transport species.
      character(len=name_len) :: species_name(max_species) = ""   !< List of species names.
      real(rk) :: species_diffusivity(max_species) = 0.0_rk       !< Constant species diffusivity \(D_k\) [m^2/s].
      real(rk) :: initial_Y(max_species) = 0.0_rk                 !< Legacy index/name-matched global initial mass fractions.
      character(len=1024) :: initial_composition = ""              !< Preferred named initial mixture, e.g. "CH4:0.028,O2:0.225,N2:0.747".

      !> @name Operator-split Cantera chemistry controls
      integer :: chemistry_update_interval = 1                    !< Chemistry source update cadence [steps].
      real(rk) :: chemistry_rtol = 1.0e-8_rk                      !< Relative tolerance for Cantera ReactorNet chemistry integration.
      real(rk) :: chemistry_atol = 1.0e-15_rk                     !< Absolute tolerance for Cantera ReactorNet chemistry integration.
      integer :: chemistry_max_steps = 10000                     !< Maximum internal ReactorNet steps per cell chemistry advance.
      logical :: chemistry_energy_enabled = .true.                !< If true, Cantera chemistry solves the reactor energy equation.
      real(rk) :: chemistry_temperature_cutoff = zero             !< Skip chemistry below this temperature [K]; zero disables the cutoff.
      real(rk) :: chemistry_min_reactive_mass_fraction = zero     !< Legacy skip if max(Y_k) is below this threshold; zero disables.
      integer :: chemistry_n_active_species = 0                  !< Number of named species used for optional chemistry activity screening.
      character(len=name_len) :: chemistry_active_species_name(max_species) = "" !< Species names used for optional chemistry screening.
      real(rk) :: chemistry_active_species_threshold = zero      !< Skip chemistry if all named active-species Y are below this threshold.
      real(rk) :: chemistry_max_dT_per_step = zero              !< Optional chemistry dt/source limiter threshold [K]; zero disables.
      real(rk) :: chemistry_max_rel_rho_change_per_step = zero  !< Optional chemistry density-change threshold; zero disables.
      real(rk) :: chemistry_max_dY_per_step = zero              !< Optional chemistry species-change threshold; zero disables.
      real(rk) :: chemistry_dt_safety = 0.5_rk                  !< Safety factor for next-step chemistry dt cuts.
      real(rk) :: chemistry_dt_min_factor = 0.1_rk              !< Minimum multiplicative dt cut in one chemistry-control event.
      logical :: chemistry_limit_source_update = .false.        !< If true, under-relax/clip committed chemistry source update.
      real(rk) :: chemistry_source_relaxation = one             !< Max fraction of Cantera chemistry update to commit when source limiting.
      logical :: enable_chemistry_load_balancing = .false.      !< If true, distribute active chemistry cells globally round-robin.

      !> @name Chemistry and Reacting Flow Stability Upgrades
      logical :: enable_chemistry_subcycling = .false.          !< If true, subcycle Cantera chemistry integration.
      integer :: max_chemistry_subcycles = 8                    !< Upper limit on chemistry subcycles per timestep.
      real(rk) :: subcycling_dT_threshold = 10.0_rk             !< Local temperature change threshold [K] to trigger subcycling.
      logical :: reject_on_nan = .true.                         !< If true, reject trial steps containing NaN/Inf values.
      real(rk) :: min_density_ratio_limit = 0.05_rk             !< Minimum density ratio relative to background before step rejection.
      real(rk) :: max_temperature_limit = 3200.0_rk             !< Maximum temperature limit [K] before step rejection.
      real(rk) :: min_temperature_limit = 200.0_rk              !< Minimum temperature limit [K] before step rejection.
      real(rk) :: max_pressure_residual_limit = 1.0e-2_rk       !< Maximum allowed PCG pressure residual before step rejection.

      !> @name Ignition-kernel initialization
      logical :: enable_ignition_kernel = .false.                 !< If true, overwrite an initial hot spherical kernel before the first step.
      character(len=name_len) :: ignition_kernel_shape = "sphere" !< First supported shape is sphere.
      real(rk) :: ignition_kernel_center(3) = zero                !< Kernel center [m].
      real(rk) :: ignition_kernel_radius = zero                   !< Kernel radius [m].
      real(rk) :: ignition_kernel_T = zero                        !< Kernel temperature [K].
      character(len=1024) :: ignition_kernel_composition = ""  !< Optional named composition inside kernel.
      character(len=name_len) :: ignition_kernel_blend = "overwrite" !< First supported blend mode is overwrite.
      
      !> @name Internal Registry (Reacting discovery)
      integer :: namelist_nspecies = 0                            !< Number of species specified in `case.nml`.
      character(len=name_len) :: namelist_species_name(max_species) = "" !< Names specified in `case.nml`.
      
      !> @name Cantera Bridge Integration
      logical :: enable_cantera_fluid = .false.                   !< Use Cantera for mixture-averaged viscosity; flow density remains `rho`.
      logical :: enable_cantera_species = .false.                 !< Use Cantera for mixture-averaged \(D_k\).
      character(len=path_len) :: cantera_mech_file = "gri30.yaml" !< Path to YAML/CTI mechanism file.
      character(len=name_len) :: cantera_phase_name = ""          !< Optional Cantera phase name; blank uses first/default phase.
      integer :: cantera_reaction_count = 0                       !< Number of reactions reported by the selected Cantera phase.
      real(rk) :: background_temp = 300.0_rk                      !< Fixed temperature for property evaluation [K].
      real(rk) :: background_press = 101325.0_rk                  !< Uniform thermodynamic pressure \(p_0\) passed to Cantera [Pa].

      !> @name Enthalpy Energy Equation Controls
      logical :: enable_energy = .false.                         !< Enable enthalpy/temperature field storage.
      logical :: enable_cantera_thermo = .false.                 !< Use Cantera for sensible h(T,Y,p0) and T(h,Y,p0).
      logical :: enable_species_enthalpy_diffusion = .false.     !< Include optional species-enthalpy diffusion correction in the \(h\) equation.
      integer :: thermo_update_interval = 1                     !< Reserved Cantera thermo interval [steps]; must remain 1.
      character(len=name_len) :: thermo_default_species = 'N2'   !< Single-species Cantera thermo fallback when species transport is off.
      real(rk) :: initial_T = 300.0_rk                           !< Initial gas temperature [K].
      real(rk) :: energy_reference_T = 298.15_rk                 !< Reference temperature for sensible enthalpy [K].
      real(rk) :: energy_reference_h = zero                      !< Reference sensible enthalpy for constant-cp mode [J/kg].
      real(rk) :: energy_cp = 1005.0_rk                          !< Constant heat capacity for non-Cantera thermo [J/kg/K].
      real(rk) :: energy_lambda = 2.6e-2_rk                      !< Constant thermal conductivity for non-Cantera thermo [W/m/K].

      !> @name Spectral radiation coupling controls
      logical :: enable_radiation = .false.                       !< Enable radiation-source updates on the radiation communicator.
      integer :: radiation_update_interval = 1                    !< Recompute qrad every N flow steps; otherwise reuse last qrad.
      integer :: radiation_n_wavenumbers = 0                      !< Number of spectral/wavenumber tasks decomposed over radiation MPI ranks.
      integer :: radiation_n_species = 0                          !< Number of transported species used by the radiation model.
      character(len=name_len) :: radiation_species_name(max_species) = "" !< Radiation-relevant transported species names.
      character(len=name_len) :: radiation_source_model = "none" !< Radiation source model: none or spectral_test.
      real(rk) :: radiation_source_scale = 0.0_rk                 !< Test-model source scale [W/m^3].
      integer :: radiation_n_scalars = 0                            !< Reserved number of generic scalar arrays for radiation.
      character(len=name_len) :: radiation_scalar_name(max_species) = "" !< Reserved selected scalar names for future radiation models.
      character(len=name_len) :: radiation_pressure_source = "background" !< Pressure source for radiation state: background or system.
      logical :: radiation_debug = .false.                          !< Write extra gather/checksum radiation diagnostics.
      logical :: write_radiation_diagnostics = .true.               !< Write radiation_setup.txt and radiation_diagnostics.csv.
      real(rk) :: radiation_absorption_coeff = 1.0_rk               !< Constant grey gas absorption coefficient [1/m].
      real(rk) :: radiation_emissivity = 1.0_rk                     !< Wall boundary surface emissivity [0-1].
      character(len=name_len) :: radiation_dom_quadrature = "s2"    !< DOM angular quadrature type: s2 or s4.


      !> @name Profiling Controls
      logical :: enable_profiling = .true.                       !< Enable wall-clock profiling.
      logical :: nested_profiling = .true.                       !< If true, print nested profiling tree.
      logical :: write_cantera_cache_stats = .false.              !< If true, print final Cantera cache hit/miss statistics.
      logical :: write_chemistry_screening_stats = .false.        !< If true, print final chemistry screening/skipped-cell statistics.
      logical :: variable_density_debug = .false.                 !< If true, print verbose experimental variable-density diagnostics.
      integer :: variable_density_debug_interval = 1              !< Print variable-density debug every N energy steps when enabled.
   end type case_params_t

   public :: read_case_params

contains

   !> Orchestrates the reading of all namelist blocks from the configuration file.
   !!
   !! Performs a sequential read of all expected namelist groups and 
   !! triggers a final validation pass to ensure physical consistency.
   !!
   !! @param filename Path to the `.nml` file (usually `case.nml`).
   !! @param params Container to be populated with parsed values.
   subroutine read_case_params(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      call read_mesh_input(filename, params)
      call read_time_input(filename, params)
      call read_fluid_input(filename, params)
      call read_solver_input(filename, params)
      call read_boundary_input(filename, params)
      call read_species_input(filename, params)
      call read_chemistry_input(filename, params)
      call read_energy_input(filename, params)
      call read_ignition_input(filename, params)
      call read_radiation_input(filename, params)
      call read_output_input(filename, params)
      call read_restart_input(filename, params)
      call read_profiling_input(filename, params)
      call validate_params(params)

      ! Density semantics:
      ! constant-density mode keeps transport%rho=params%rho and treats
      ! energy%rho_thermo as a diagnostic; variable-density density_eos="cantera"
      ! requires Cantera thermo and syncs transport%rho from energy%rho_thermo.
      ! The selected YAML phase defines the actual EOS.  "ideal_gas" is parsed
      ! only to reserve the name and is rejected for active variable-density use.
      params%density_eos = trim(lowercase(params%density_eos))
      if (len_trim(params%density_eos) == 0) params%density_eos = "constant"

      select case (trim(params%density_eos))
      case ("constant", "cantera", "ideal_gas")
         continue
      case default
         call fatal_error("input", "density_eos must be one of: constant, cantera, ideal_gas")
      end select

      if (params%enable_variable_density) then
         select case (trim(params%density_eos))
         case ("cantera")
            if (.not. params%enable_energy) then
               call fatal_error("input", "enable_variable_density with density_eos=cantera requires enable_energy=.true.")
            end if
            if (.not. params%enable_cantera_thermo) then
               call fatal_error("input", "enable_variable_density with density_eos=cantera requires enable_cantera_thermo=.true.")
            end if
            if (params%thermo_update_interval /= 1) then
               call fatal_error("input", "enable_variable_density requires thermo_update_interval=1")
            end if
         case ("ideal_gas")
            call fatal_error("input", "density_eos=ideal_gas is parsed but not implemented for variable-density mode yet")
         case default
            call fatal_error("input", "enable_variable_density currently requires density_eos=cantera")
         end select
      end if

   end subroutine read_case_params


   !> Reads the `&mesh_input` namelist block.
   subroutine read_mesh_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      character(len=path_len) :: mesh_dir
      integer :: unit_id, ios

      namelist /mesh_input/ mesh_dir

      mesh_dir = params%mesh_dir

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=mesh_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &mesh_input')

      params%mesh_dir = mesh_dir
   end subroutine read_mesh_input


   !> Reads the `&time_input` namelist block.
   subroutine read_time_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      integer :: nsteps, output_interval
      real(rk) :: dt, max_cfl, max_dt, min_dt, dt_growth_limit
      logical :: use_dynamic_dt
      logical :: enable_step_rejection, reject_on_pressure_max_iter, stop_on_min_dt_failure
      integer :: max_step_retries
      real(rk) :: step_reject_cut_factor, step_reject_min_dt
      real(rk) :: max_KE_growth_per_step, ke_reject_floor, max_cfl_overshoot_factor
      integer :: unit_id, ios

      namelist /time_input/ nsteps, dt, output_interval, use_dynamic_dt, max_cfl, &
                             max_dt, min_dt, dt_growth_limit, &
                             enable_step_rejection, max_step_retries, step_reject_cut_factor, &
                             step_reject_min_dt, max_KE_growth_per_step, ke_reject_floor, max_cfl_overshoot_factor, &
                             reject_on_pressure_max_iter, stop_on_min_dt_failure

      nsteps = params%nsteps
      dt = params%dt
      output_interval = params%output_interval
      use_dynamic_dt = params%use_dynamic_dt
      max_cfl = params%max_cfl
      max_dt = params%max_dt
      min_dt = params%min_dt
      dt_growth_limit = params%dt_growth_limit
      enable_step_rejection = params%enable_step_rejection
      max_step_retries = params%max_step_retries
      step_reject_cut_factor = params%step_reject_cut_factor
      step_reject_min_dt = params%step_reject_min_dt
      max_KE_growth_per_step = params%max_KE_growth_per_step
      ke_reject_floor = params%ke_reject_floor
      max_cfl_overshoot_factor = params%max_cfl_overshoot_factor
      reject_on_pressure_max_iter = params%reject_on_pressure_max_iter
      stop_on_min_dt_failure = params%stop_on_min_dt_failure

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=time_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &time_input')

      params%nsteps = nsteps
      params%dt = dt
      params%output_interval = output_interval
      params%use_dynamic_dt = use_dynamic_dt
      params%max_cfl = max_cfl
      params%max_dt = max_dt
      params%min_dt = min_dt
      params%dt_growth_limit = dt_growth_limit
      params%enable_step_rejection = enable_step_rejection
      params%max_step_retries = max_step_retries
      params%step_reject_cut_factor = step_reject_cut_factor
      params%step_reject_min_dt = step_reject_min_dt
      params%max_KE_growth_per_step = max_KE_growth_per_step
      params%ke_reject_floor = ke_reject_floor
      params%max_cfl_overshoot_factor = max_cfl_overshoot_factor
      params%reject_on_pressure_max_iter = reject_on_pressure_max_iter
      params%stop_on_min_dt_failure = stop_on_min_dt_failure
   end subroutine read_time_input


   !> Reads the `&fluid_input` namelist block.
   subroutine read_fluid_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      real(rk) :: rho, nu
      logical :: enable_cantera
      character(len=path_len) :: cantera_mech_file
      character(len=name_len) :: cantera_phase_name
      real(rk) :: background_temp
      real(rk) :: background_press
      integer :: transport_update_interval
      integer :: unit_id, ios
      logical :: enable_variable_density
      logical :: enable_density_corrector_projection
      logical :: upwind_dilatation_density
      character(len=name_len) :: density_eos = "constant" !< Density EOS selector: constant, cantera, ideal_gas.
      logical :: enable_variable_nu
      namelist /fluid_input/ rho, nu, enable_cantera, cantera_mech_file, cantera_phase_name, &
                              background_temp, background_press, transport_update_interval, &
                              enable_variable_density, enable_density_corrector_projection, &
                              density_eos, enable_variable_nu, upwind_dilatation_density

      rho = params%rho
      nu = params%nu
      enable_cantera = params%enable_cantera_fluid
      cantera_mech_file = params%cantera_mech_file
      cantera_phase_name = params%cantera_phase_name
      background_temp = params%background_temp
      background_press = params%background_press
      transport_update_interval = params%transport_update_interval
      enable_variable_density = params%enable_variable_density
      enable_density_corrector_projection = params%enable_density_corrector_projection
      density_eos = params%density_eos
      enable_variable_nu = params%enable_variable_nu
      upwind_dilatation_density = params%upwind_dilatation_density

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=fluid_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &fluid_input')

      params%rho = rho
      params%nu = nu
      params%enable_cantera_fluid = enable_cantera
      params%enable_variable_density = enable_variable_density
      params%enable_density_corrector_projection = enable_density_corrector_projection
      params%density_eos = trim(density_eos)
      params%enable_variable_nu = enable_variable_nu
      params%upwind_dilatation_density = upwind_dilatation_density
      params%cantera_mech_file = cantera_mech_file
      params%cantera_phase_name = trim(cantera_phase_name)
      params%background_temp = background_temp
      params%background_press = background_press
      params%transport_update_interval = transport_update_interval
   end subroutine read_fluid_input


   !> Reads the `&solver_input` namelist block.
   subroutine read_solver_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      integer :: pressure_max_iter
      real(rk) :: pressure_tol
      real(rk) :: pressure_abs_tol
      real(rk) :: body_force_x, body_force_y, body_force_z
      character(len=name_len) :: convection_scheme
      character(len=name_len) :: momentum_convection_scheme
      character(len=name_len) :: species_convection_scheme
      character(len=name_len) :: energy_convection_scheme
      character(len=name_len) :: scalar_limiter
      character(len=name_len) :: species_time_scheme
      character(len=name_len) :: energy_time_scheme
      integer :: unit_id, ios

      namelist /solver_input/ pressure_max_iter, pressure_tol, pressure_abs_tol, &
                              body_force_x, body_force_y, body_force_z, &
                              convection_scheme, momentum_convection_scheme, &
                              species_convection_scheme, energy_convection_scheme, &
                              scalar_limiter, species_time_scheme, energy_time_scheme

      pressure_max_iter = params%pressure_max_iter
      pressure_tol = params%pressure_tol
      pressure_abs_tol = params%pressure_abs_tol
      body_force_x = params%body_force(1)
      body_force_y = params%body_force(2)
      body_force_z = params%body_force(3)
      convection_scheme = params%convection_scheme
      momentum_convection_scheme = params%momentum_convection_scheme
      species_convection_scheme = params%species_convection_scheme
      energy_convection_scheme = params%energy_convection_scheme
      scalar_limiter = params%scalar_limiter
      species_time_scheme = params%species_time_scheme
      energy_time_scheme = params%energy_time_scheme

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=solver_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &solver_input')

      params%pressure_max_iter = pressure_max_iter
      params%pressure_tol = pressure_tol
      params%pressure_abs_tol = pressure_abs_tol
      params%body_force = [body_force_x, body_force_y, body_force_z]
      params%convection_scheme = trim(lowercase(convection_scheme))
      params%momentum_convection_scheme = trim(lowercase(momentum_convection_scheme))
      if (len_trim(params%momentum_convection_scheme) == 0) then
         params%momentum_convection_scheme = params%convection_scheme
      end if
      ! The legacy field remains the value consumed by the existing projection solver.
      params%convection_scheme = params%momentum_convection_scheme
      params%species_convection_scheme = trim(lowercase(species_convection_scheme))
      params%energy_convection_scheme = trim(lowercase(energy_convection_scheme))
      params%scalar_limiter = trim(lowercase(scalar_limiter))
      params%species_time_scheme = trim(lowercase(species_time_scheme))
      params%energy_time_scheme = trim(lowercase(energy_time_scheme))
   end subroutine read_solver_input


   !> Reads the `&boundary_input` namelist block.
   subroutine read_boundary_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      integer :: n_patches
      character(len=name_len) :: patch_name(max_patches)
      character(len=name_len) :: patch_type(max_patches)
      character(len=name_len) :: patch_velocity_type(max_patches)
      character(len=name_len) :: patch_pressure_type(max_patches)
      character(len=name_len) :: patch_temperature_type(max_patches)
      character(len=name_len) :: patch_species_type(max_patches)
      real(rk) :: patch_u(max_patches)
      real(rk) :: patch_v(max_patches)
      real(rk) :: patch_w(max_patches)
      real(rk) :: patch_mass_flux(max_patches)
      real(rk) :: patch_p(max_patches)
      real(rk) :: patch_dpdn(max_patches)
      real(rk) :: patch_T(max_patches)
      real(rk), save :: patch_Y(max_species, max_patches) = zero
      character(len=1024), save :: patch_composition(max_patches) = ""
      integer :: unit_id, ios, i

      namelist /boundary_input/ n_patches, patch_name, patch_type, &
                                 patch_velocity_type, patch_pressure_type, &
                                 patch_temperature_type, patch_species_type, &
                                 patch_u, patch_v, patch_w, patch_mass_flux, &
                                 patch_p, patch_dpdn, patch_T, patch_Y, patch_composition

      n_patches = params%n_patches
      patch_name = params%patch_name
      patch_type = params%patch_type
      patch_velocity_type = params%patch_velocity_type
      patch_pressure_type = params%patch_pressure_type
      patch_temperature_type = params%patch_temperature_type
      patch_species_type = params%patch_species_type
      patch_u = params%patch_u
      patch_v = params%patch_v
      patch_w = params%patch_w
      patch_mass_flux = params%patch_mass_flux
      patch_p = params%patch_p
      patch_dpdn = params%patch_dpdn
      patch_T = params%patch_T
      patch_Y = params%patch_Y
      patch_composition = params%patch_composition

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=boundary_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &boundary_input')

      params%n_patches = n_patches
      params%patch_name = patch_name
      params%patch_type = patch_type
      params%patch_velocity_type = patch_velocity_type
      params%patch_pressure_type = patch_pressure_type
      params%patch_temperature_type = patch_temperature_type
      params%patch_species_type = patch_species_type

      do i = 1, max_patches
         params%patch_type(i) = trim(lowercase(params%patch_type(i)))
         params%patch_velocity_type(i) = trim(lowercase(params%patch_velocity_type(i)))
         params%patch_pressure_type(i) = trim(lowercase(params%patch_pressure_type(i)))
         params%patch_temperature_type(i) = trim(lowercase(params%patch_temperature_type(i)))
         params%patch_species_type(i) = trim(lowercase(params%patch_species_type(i)))
      end do

      params%patch_u = patch_u
      params%patch_v = patch_v
      params%patch_w = patch_w
      params%patch_mass_flux = patch_mass_flux
      params%patch_p = patch_p
      params%patch_dpdn = patch_dpdn
      params%patch_T = patch_T
      params%patch_Y = patch_Y
      params%patch_composition = patch_composition
   end subroutine read_boundary_input


   !> Reads the `&profiling_input` block.
   !!
   !! This parser is intentionally explicit because profiling is optional and
   !! should be robust to block ordering in case.nml.
   subroutine read_profiling_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      character(len=512) :: line
      character(len=256) :: key, value
      integer :: unit_id, ios, eqpos, comment_pos, comma_pos
      logical :: in_block, found_block

      in_block = .false.
      found_block = .false.

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)
      if (ios /= 0) return

      do
         read(unit_id, '(a)', iostat=ios) line
         if (ios /= 0) exit

         comment_pos = index(line, '!')
         if (comment_pos > 0) line = line(:comment_pos-1)

         line = adjustl(line)
         if (len_trim(line) == 0) cycle

         if (.not. in_block) then
            if (index(lowercase(line), '&profiling_input') == 1) then
               in_block = .true.
               found_block = .true.
            end if
            cycle
         end if

         if (index(line, '/') > 0) exit

         eqpos = index(line, '=')
         if (eqpos <= 0) cycle

         key = trim(adjustl(lowercase(line(:eqpos-1))))
         value = trim(adjustl(lowercase(line(eqpos+1:))))

         comma_pos = index(value, ',')
         if (comma_pos > 0) value = value(:comma_pos-1)
         value = trim(value)

         select case (trim(key))
         case ('enable_profiling')
            call parse_logical_value(value, params%enable_profiling, 'enable_profiling')
         case ('nested_profiling')
            call parse_logical_value(value, params%nested_profiling, 'nested_profiling')
         case ('write_cantera_cache_stats')
            call parse_logical_value(value, params%write_cantera_cache_stats, 'write_cantera_cache_stats')
         case ('write_chemistry_screening_stats')
            call parse_logical_value(value, params%write_chemistry_screening_stats, 'write_chemistry_screening_stats')
         case ('variable_density_debug')
            call parse_logical_value(value, params%variable_density_debug, 'variable_density_debug')
         case ('variable_density_debug_interval')
            read(value, *, iostat=ios) params%variable_density_debug_interval
            if (ios /= 0) call fatal_error('input', 'invalid integer for variable_density_debug_interval: '//trim(value))
         case default
            call fatal_error('input', 'unknown variable in &profiling_input: '//trim(key))
         end select
      end do

      close(unit_id)

      if (found_block .and. in_block .and. ios /= 0) then
         call fatal_error('input', 'unterminated &profiling_input block')
      end if

   contains

      !> Parses a string representation of a logical value into a Fortran logical.
      !!
      !! Supports common formats like '.true.', 'true', 't', and their false counterparts.
      !!
      !! @param value_text The string to parse.
      !! @param output_value The resulting logical value.
      !! @param field_name The name of the field being parsed (used for error reporting).
      subroutine parse_logical_value(value_text, output_value, field_name)
         character(len=*), intent(in) :: value_text
         logical, intent(out) :: output_value
         character(len=*), intent(in) :: field_name

         select case (trim(value_text))
         case ('.true.', 'true', 't', '.t.')
            output_value = .true.
         case ('.false.', 'false', 'f', '.f.')
            output_value = .false.
         case default
            call fatal_error('input', 'invalid logical for '//trim(field_name)//': '//trim(value_text))
         end select
      end subroutine parse_logical_value

   end subroutine read_profiling_input



   !> Reads the optional `&chemistry_input` namelist block.
   !!
   !! Chemistry is operator-split from advection/diffusion and is solved by the
   !! Cantera C++ bridge as a cell-local constant-pressure reactor.  The block
   !! is optional; defaults are conservative and do not enable chemistry by
   !! themselves.  `enable_reactions` remains in `&species_input` for backward
   !! compatibility with existing cases.
   subroutine read_chemistry_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      integer :: chemistry_update_interval
      real(rk) :: chemistry_rtol, chemistry_atol
      integer :: chemistry_max_steps
      logical :: chemistry_energy_enabled
      real(rk) :: chemistry_temperature_cutoff
      real(rk) :: chemistry_min_reactive_mass_fraction
      integer :: chemistry_n_active_species
      character(len=name_len) :: chemistry_active_species_name(max_species)
      real(rk) :: chemistry_active_species_threshold
      real(rk) :: chemistry_max_dT_per_step, chemistry_max_rel_rho_change_per_step
      real(rk) :: chemistry_max_dY_per_step, chemistry_dt_safety, chemistry_dt_min_factor
      logical :: chemistry_limit_source_update
      real(rk) :: chemistry_source_relaxation
      logical :: enable_chemistry_subcycling
      integer :: max_chemistry_subcycles
      real(rk) :: subcycling_dT_threshold
      logical :: reject_on_nan
      real(rk) :: min_density_ratio_limit
      real(rk) :: max_temperature_limit
      real(rk) :: min_temperature_limit
      real(rk) :: max_pressure_residual_limit
      logical :: enable_chemistry_load_balancing
      integer :: unit_id, ios, k

      namelist /chemistry_input/ chemistry_update_interval, chemistry_rtol, chemistry_atol, &
                                  chemistry_max_steps, chemistry_energy_enabled, &
                                  chemistry_temperature_cutoff, chemistry_min_reactive_mass_fraction, &
                                  chemistry_n_active_species, chemistry_active_species_name, &
                                  chemistry_active_species_threshold, &
                                  chemistry_max_dT_per_step, chemistry_max_rel_rho_change_per_step, &
                                  chemistry_max_dY_per_step, chemistry_dt_safety, chemistry_dt_min_factor, &
                                  chemistry_limit_source_update, chemistry_source_relaxation, &
                                  enable_chemistry_subcycling, max_chemistry_subcycles, subcycling_dT_threshold, &
                                  reject_on_nan, min_density_ratio_limit, max_temperature_limit, &
                                  min_temperature_limit, max_pressure_residual_limit, &
                                  enable_chemistry_load_balancing

      chemistry_update_interval = params%chemistry_update_interval
      chemistry_rtol = params%chemistry_rtol
      chemistry_atol = params%chemistry_atol
      chemistry_max_steps = params%chemistry_max_steps
      chemistry_energy_enabled = params%chemistry_energy_enabled
      chemistry_temperature_cutoff = params%chemistry_temperature_cutoff
      chemistry_min_reactive_mass_fraction = params%chemistry_min_reactive_mass_fraction
      chemistry_n_active_species = params%chemistry_n_active_species
      chemistry_active_species_name = params%chemistry_active_species_name
      chemistry_active_species_threshold = params%chemistry_active_species_threshold
      chemistry_max_dT_per_step = params%chemistry_max_dT_per_step
      chemistry_max_rel_rho_change_per_step = params%chemistry_max_rel_rho_change_per_step
      chemistry_max_dY_per_step = params%chemistry_max_dY_per_step
      chemistry_dt_safety = params%chemistry_dt_safety
      chemistry_dt_min_factor = params%chemistry_dt_min_factor
      chemistry_limit_source_update = params%chemistry_limit_source_update
      chemistry_source_relaxation = params%chemistry_source_relaxation
      enable_chemistry_subcycling = params%enable_chemistry_subcycling
      max_chemistry_subcycles = params%max_chemistry_subcycles
      subcycling_dT_threshold = params%subcycling_dT_threshold
      reject_on_nan = params%reject_on_nan
      min_density_ratio_limit = params%min_density_ratio_limit
      max_temperature_limit = params%max_temperature_limit
      min_temperature_limit = params%min_temperature_limit
      max_pressure_residual_limit = params%max_pressure_residual_limit
      enable_chemistry_load_balancing = params%enable_chemistry_load_balancing

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=chemistry_input, iostat=ios)
         close(unit_id)
      end if

      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &chemistry_input. Check for unknown variables or typos.')
      end if

      if (ios == 0) then
         params%chemistry_update_interval = chemistry_update_interval
         params%chemistry_rtol = chemistry_rtol
         params%chemistry_atol = chemistry_atol
         params%chemistry_max_steps = chemistry_max_steps
         params%chemistry_energy_enabled = chemistry_energy_enabled
         params%chemistry_temperature_cutoff = chemistry_temperature_cutoff
         params%chemistry_min_reactive_mass_fraction = chemistry_min_reactive_mass_fraction
         params%chemistry_n_active_species = chemistry_n_active_species
         params%chemistry_active_species_name = chemistry_active_species_name
         params%chemistry_active_species_threshold = chemistry_active_species_threshold
         params%chemistry_max_dT_per_step = chemistry_max_dT_per_step
         params%chemistry_max_rel_rho_change_per_step = chemistry_max_rel_rho_change_per_step
         params%chemistry_max_dY_per_step = chemistry_max_dY_per_step
         params%chemistry_dt_safety = chemistry_dt_safety
         params%chemistry_dt_min_factor = chemistry_dt_min_factor
         params%chemistry_limit_source_update = chemistry_limit_source_update
         params%chemistry_source_relaxation = chemistry_source_relaxation
         params%enable_chemistry_subcycling = enable_chemistry_subcycling
         params%max_chemistry_subcycles = max_chemistry_subcycles
         params%subcycling_dT_threshold = subcycling_dT_threshold
         params%reject_on_nan = reject_on_nan
         params%min_density_ratio_limit = min_density_ratio_limit
         params%max_temperature_limit = max_temperature_limit
         params%min_temperature_limit = min_temperature_limit
         params%max_pressure_residual_limit = max_pressure_residual_limit
         params%enable_chemistry_load_balancing = enable_chemistry_load_balancing

         if (params%chemistry_n_active_species <= 0) then
            do k = 1, max_species
               if (len_trim(params%chemistry_active_species_name(k)) > 0) then
                  params%chemistry_n_active_species = k
               end if
            end do
         end if
      end if
   end subroutine read_chemistry_input


   !> Reads the optional `&radiation_input` namelist block.
   !!
   !! Radiation is intentionally configured independently from flow MPI.  Flow
   !! remains spatially decomposed while radiation decomposes spectral tasks
   !! (`radiation_n_wavenumbers`) over the radiation communicator and evaluates
   !! source terms on the full replicated mesh.
   subroutine read_radiation_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      logical :: enable_radiation
      integer :: radiation_update_interval, radiation_n_wavenumbers, radiation_n_species
      integer :: radiation_n_scalars
      character(len=name_len), save :: radiation_species_name(max_species) = ""
      character(len=name_len), save :: radiation_scalar_name(max_species) = ""
      character(len=name_len) :: radiation_source_model, radiation_pressure_source, radiation_dom_quadrature
      real(rk) :: radiation_source_scale
      logical :: radiation_debug, write_radiation_diagnostics
      real(rk) :: radiation_absorption_coeff, radiation_emissivity
      integer :: unit_id, ios, i
 
      namelist /radiation_input/ enable_radiation, radiation_update_interval, &
                                  radiation_n_wavenumbers, radiation_n_species, &
                                  radiation_species_name, radiation_source_model, &
                                  radiation_source_scale, radiation_n_scalars, &
                                  radiation_scalar_name, radiation_pressure_source, &
                                  radiation_debug, write_radiation_diagnostics, &
                                  radiation_absorption_coeff, radiation_emissivity, &
                                  radiation_dom_quadrature
 
      enable_radiation = params%enable_radiation
      radiation_update_interval = params%radiation_update_interval
      radiation_n_wavenumbers = params%radiation_n_wavenumbers
      radiation_n_species = params%radiation_n_species
      radiation_species_name = params%radiation_species_name
      radiation_scalar_name = params%radiation_scalar_name
      radiation_source_model = params%radiation_source_model
      radiation_source_scale = params%radiation_source_scale
      radiation_n_scalars = params%radiation_n_scalars
      radiation_pressure_source = params%radiation_pressure_source
      radiation_debug = params%radiation_debug
      write_radiation_diagnostics = params%write_radiation_diagnostics
      radiation_absorption_coeff = params%radiation_absorption_coeff
      radiation_emissivity = params%radiation_emissivity
      radiation_dom_quadrature = params%radiation_dom_quadrature
 
      call open_namelist_file(filename, unit_id, ios)
 
      if (ios == 0) then
         read(unit_id, nml=radiation_input, iostat=ios)
         close(unit_id)
      end if
 
      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &radiation_input. Check for unknown variables or typos.')
      end if
 
      if (ios == 0) then
         params%enable_radiation = enable_radiation
         params%radiation_update_interval = radiation_update_interval
         params%radiation_n_wavenumbers = radiation_n_wavenumbers
         params%radiation_n_species = radiation_n_species
         do i = 1, max_species
            params%radiation_species_name(i) = trim(radiation_species_name(i))
         end do
         params%radiation_source_model = trim(lowercase(radiation_source_model))
         params%radiation_source_scale = radiation_source_scale
         params%radiation_n_scalars = radiation_n_scalars
         do i = 1, max_species
            params%radiation_scalar_name(i) = trim(radiation_scalar_name(i))
         end do
         params%radiation_pressure_source = trim(lowercase(radiation_pressure_source))
         params%radiation_debug = radiation_debug
         params%write_radiation_diagnostics = write_radiation_diagnostics
         params%radiation_absorption_coeff = radiation_absorption_coeff
         params%radiation_emissivity = radiation_emissivity
         params%radiation_dom_quadrature = trim(lowercase(radiation_dom_quadrature))
      end if
   end subroutine read_radiation_input


   !> Reads the `&output_input` namelist block.
   !> Reads the `&energy_input` namelist block.
   !!
   !! Reads controls for sensible-enthalpy transport. Cantera thermo, when
   !! enabled, preserves transported `h` and recovers `T(h,Y,p0)` every energy
   !! step; `thermo_update_interval` is reserved and must remain 1.
   subroutine read_energy_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      logical :: enable_energy
      logical :: enable_cantera_thermo
      logical :: enable_species_enthalpy_diffusion
      integer :: thermo_update_interval
      character(len=name_len) :: thermo_default_species
      real(rk) :: initial_T
      real(rk) :: energy_reference_T
      real(rk) :: energy_reference_h
      real(rk) :: energy_cp
      real(rk) :: energy_lambda
      integer :: unit_id, ios

      namelist /energy_input/ enable_energy, enable_cantera_thermo, enable_species_enthalpy_diffusion, &
                               thermo_default_species, thermo_update_interval, &
                               initial_T, energy_reference_T, energy_reference_h, &
                               energy_cp, energy_lambda

      enable_energy = params%enable_energy
      enable_cantera_thermo = params%enable_cantera_thermo
      enable_species_enthalpy_diffusion = params%enable_species_enthalpy_diffusion
      thermo_update_interval = params%thermo_update_interval
      thermo_default_species = params%thermo_default_species
      initial_T = params%initial_T
      energy_reference_T = params%energy_reference_T
      energy_reference_h = params%energy_reference_h
      energy_cp = params%energy_cp
      energy_lambda = params%energy_lambda

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=energy_input, iostat=ios)
         close(unit_id)
      end if

      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &energy_input. Check for unknown variables or typos.')
      end if

      if (ios == 0) then
         params%enable_energy = enable_energy
         params%enable_cantera_thermo = enable_cantera_thermo
         params%enable_species_enthalpy_diffusion = enable_species_enthalpy_diffusion
         params%thermo_update_interval = thermo_update_interval
         params%thermo_default_species = trim(thermo_default_species)
         params%initial_T = initial_T
         params%energy_reference_T = energy_reference_T
         params%energy_reference_h = energy_reference_h
         params%energy_cp = energy_cp
         params%energy_lambda = energy_lambda
      end if
   end subroutine read_energy_input


   !> Reads the optional `&ignition_input` namelist block.
   !!
   !! The ignition kernel is an initialization aid, not a time-dependent source.
   !! It overwrites temperature and optionally composition inside a spherical
   !! region before the first timestep, then the caller refreshes thermodynamics.
   subroutine read_ignition_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      logical :: enable_ignition_kernel
      character(len=name_len) :: ignition_kernel_shape
      real(rk) :: ignition_kernel_center(3)
      real(rk) :: ignition_kernel_radius
      real(rk) :: ignition_kernel_T
      character(len=1024) :: ignition_kernel_composition
      character(len=name_len) :: ignition_kernel_blend
      integer :: unit_id, ios

      namelist /ignition_input/ enable_ignition_kernel, ignition_kernel_shape, &
                                 ignition_kernel_center, ignition_kernel_radius, &
                                 ignition_kernel_T, ignition_kernel_composition, &
                                 ignition_kernel_blend

      enable_ignition_kernel = params%enable_ignition_kernel
      ignition_kernel_shape = params%ignition_kernel_shape
      ignition_kernel_center = params%ignition_kernel_center
      ignition_kernel_radius = params%ignition_kernel_radius
      ignition_kernel_T = params%ignition_kernel_T
      ignition_kernel_composition = params%ignition_kernel_composition
      ignition_kernel_blend = params%ignition_kernel_blend

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=ignition_input, iostat=ios)
         close(unit_id)
      end if

      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &ignition_input. Check for unknown variables or typos.')
      end if

      if (ios == 0) then
         params%enable_ignition_kernel = enable_ignition_kernel
         params%ignition_kernel_shape = trim(lowercase(ignition_kernel_shape))
         params%ignition_kernel_center = ignition_kernel_center
         params%ignition_kernel_radius = ignition_kernel_radius
         params%ignition_kernel_T = ignition_kernel_T
         params%ignition_kernel_composition = trim(ignition_kernel_composition)
         params%ignition_kernel_blend = trim(lowercase(ignition_kernel_blend))
      end if
   end subroutine read_ignition_input


   !> Reads the optional `&restart_input` block.
   subroutine read_restart_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      logical :: restart_from_file
      character(len=path_len) :: restart_file
      logical :: write_restart
      integer :: restart_interval
      character(len=path_len) :: restart_output_dir
      integer :: restart_keep
      logical :: restart_use_file_dt, restart_nsteps_are_additional
      integer :: unit_id, ios

      namelist /restart_input/ restart_from_file, restart_file, write_restart, &
                                restart_interval, restart_output_dir, restart_keep, &
                                restart_use_file_dt, restart_nsteps_are_additional

      restart_from_file = params%restart_from_file
      restart_file = params%restart_file
      write_restart = params%write_restart
      restart_interval = params%restart_interval
      restart_output_dir = params%restart_output_dir
      restart_keep = params%restart_keep
      restart_use_file_dt = params%restart_use_file_dt
      restart_nsteps_are_additional = params%restart_nsteps_are_additional

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=restart_input, iostat=ios)
         close(unit_id)
      end if

      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &restart_input. Check for unknown variables or typos.')
      end if

      if (ios == 0) then
         params%restart_from_file = restart_from_file
         params%restart_file = trim(restart_file)
         params%write_restart = write_restart
         params%restart_interval = restart_interval
         params%restart_output_dir = trim(restart_output_dir)
         params%restart_keep = restart_keep
         params%restart_use_file_dt = restart_use_file_dt
         params%restart_nsteps_are_additional = restart_nsteps_are_additional
      end if
   end subroutine read_restart_input


   subroutine read_output_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      character(len=path_len) :: output_dir
      character(len=name_len) :: terminal_detail
      logical :: write_vtu, write_diagnostics
      integer :: terminal_interval
      character(len=name_len) :: vtu_format
      integer :: unit_id, ios

      namelist /output_input/ output_dir, write_vtu, write_diagnostics, terminal_interval, terminal_detail, vtu_format

      output_dir = params%output_dir
      write_vtu = params%write_vtu
      write_diagnostics = params%write_diagnostics
      terminal_interval = params%terminal_interval
      terminal_detail = params%terminal_detail
      vtu_format = params%vtu_format

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=output_input, iostat=ios)
         close(unit_id)
      end if

      if (ios > 0) call fatal_error('input', 'failed reading &output_input')

      params%output_dir = output_dir
      params%write_vtu = write_vtu
      params%write_diagnostics = write_diagnostics
      params%terminal_interval = terminal_interval
      params%terminal_detail = trim(lowercase(terminal_detail))
      params%vtu_format = trim(lowercase(vtu_format))
   end subroutine read_output_input


   !> Reads the `&species_input` namelist block.
   !!
   !! Implements strict error checking to catch typos in chemical 
   !! species names or advection settings.
   subroutine read_species_input(filename, params)
      character(len=*), intent(in) :: filename
      type(case_params_t), intent(inout) :: params

      logical :: enable_species, enable_reactions, enable_cantera
      integer :: nspecies
      character(len=name_len), save :: species_name(max_species) = ""
      real(rk), save :: species_diffusivity(max_species) = zero
      real(rk), save :: initial_Y(max_species) = zero
      character(len=1024) :: initial_composition
      integer :: unit_id, ios

      namelist /species_input/ enable_species, enable_reactions, enable_cantera, &
                                nspecies, species_name, species_diffusivity, initial_Y, &
                                initial_composition

      enable_species = params%enable_species
      enable_reactions = params%enable_reactions
      enable_cantera = params%enable_cantera_species
      nspecies = params%nspecies
      species_name = params%species_name
      species_diffusivity = params%species_diffusivity
      initial_Y = params%initial_Y
      initial_composition = params%initial_composition

      call open_namelist_file(filename, unit_id, ios)

      if (ios == 0) then
         read(unit_id, nml=species_input, iostat=ios)
         close(unit_id)
      end if

      if (ios /= 0 .and. ios /= -1) then
         call fatal_error('input', 'failed reading &species_input. Check for unknown variables or typos.')
      end if

      if (ios == 0) then
         params%enable_species = enable_species
         params%enable_reactions = enable_reactions
         params%enable_cantera_species = enable_cantera
         params%nspecies = nspecies
         params%species_name = species_name
         params%species_diffusivity = species_diffusivity
         params%initial_Y = initial_Y
         params%initial_composition = trim(initial_composition)
         
         ! Save original namelist values
         params%namelist_nspecies = nspecies
         params%namelist_species_name = species_name
      end if
   end subroutine read_species_input


   !> Helper routine to safely open a namelist file for reading.
   subroutine open_namelist_file(filename, unit_id, ios)
      character(len=*), intent(in) :: filename
      integer, intent(out) :: unit_id
      integer, intent(out) :: ios

      open(newunit=unit_id, file=trim(filename), status='old', action='read', iostat=ios)

      if (ios /= 0) unit_id = -1
   end subroutine open_namelist_file


   !> Validates all parsed parameters against physical and algorithmic limits.
   !!
   !! Triggers `fatal_error` for non-physical values (e.g., negative density, 
   !! zero timestep) or unsupported configurations.
   subroutine validate_params(params)
      type(case_params_t), intent(in) :: params

      if (params%nsteps < 0) call fatal_error('input', 'nsteps must be non-negative')
      if (params%dt <= zero) call fatal_error('input', 'dt must be positive')
      if (params%max_dt <= zero) call fatal_error('input', 'max_dt must be positive')
      if (params%min_dt <= zero) call fatal_error('input', 'min_dt must be positive')
      if (params%min_dt > params%max_dt) call fatal_error('input', 'min_dt must be <= max_dt')
      if (params%dt_growth_limit < one) call fatal_error('input', 'dt_growth_limit must be >= 1')
      if (params%max_step_retries < 0) call fatal_error('input', 'max_step_retries must be non-negative')
      if (params%step_reject_cut_factor <= zero .or. params%step_reject_cut_factor >= one) then
         call fatal_error('input', 'step_reject_cut_factor must be in (0,1)')
      end if
      if (params%step_reject_min_dt <= zero) call fatal_error('input', 'step_reject_min_dt must be positive')
      if (params%max_KE_growth_per_step < zero) call fatal_error('input', 'max_KE_growth_per_step must be non-negative')
      if (params%ke_reject_floor < zero) call fatal_error('input', 'ke_reject_floor must be non-negative')
      if (params%max_cfl_overshoot_factor < zero) call fatal_error('input', 'max_cfl_overshoot_factor must be non-negative')
      if (params%output_interval <= 0) call fatal_error('input', 'output_interval must be positive')
      if (params%terminal_interval <= 0) call fatal_error('input', 'terminal_interval must be positive')
      select case (trim(params%terminal_detail))
      case ('brief', 'health', 'full')
         continue
      case default
         call fatal_error('input', 'terminal_detail must be one of: brief, health, full')
      end select
      if (params%rho <= zero) call fatal_error('input', 'rho must be positive')
      if (params%nu < zero) call fatal_error('input', 'nu must be non-negative')
      if (params%enable_variable_nu .and. .not. params%enable_cantera_fluid) then
         call fatal_error('input', 'enable_variable_nu=.true. requires fluid_input enable_cantera=.true.; otherwise use constant nu')
      end if
      if (params%transport_update_interval <= 0) call fatal_error('input', 'transport_update_interval must be positive')
      if (params%pressure_max_iter <= 0) call fatal_error('input', 'pressure_max_iter must be positive')
      call validate_advection_scheme_local(params%momentum_convection_scheme, 'momentum_convection_scheme', .false.)
      call validate_advection_scheme_local(params%convection_scheme, 'convection_scheme', .false.)
      call validate_advection_scheme_local(params%species_convection_scheme, 'species_convection_scheme', .true.)
      call validate_advection_scheme_local(params%energy_convection_scheme, 'energy_convection_scheme', .true.)
      call validate_limiter_name_local(params%scalar_limiter)
      call validate_time_scheme_local(params%species_time_scheme, 'species_time_scheme')
      call validate_time_scheme_local(params%energy_time_scheme, 'energy_time_scheme')

      if (params%initial_T <= zero) call fatal_error('input', 'initial_T must be positive')
      if (params%energy_reference_T <= zero) call fatal_error('input', 'energy_reference_T must be positive')
      if (params%energy_cp <= zero) call fatal_error('input', 'energy_cp must be positive')
      if (params%energy_lambda < zero) call fatal_error('input', 'energy_lambda must be non-negative')
      if (params%thermo_update_interval <= 0) call fatal_error('input', 'thermo_update_interval must be positive')
      if (params%enable_cantera_thermo .and. params%thermo_update_interval /= 1) then
         call fatal_error('input', 'thermo_update_interval values other than 1 are reserved; Cantera thermo is currently updated every energy step')
      end if
      if (params%enable_cantera_thermo .and. len_trim(params%thermo_default_species) == 0) &
         call fatal_error('input', 'thermo_default_species must not be empty when enable_cantera_thermo is true')

      if (params%enable_species_enthalpy_diffusion) then
         if (.not. params%enable_energy) then
            call fatal_error('input', 'enable_species_enthalpy_diffusion requires enable_energy=.true.')
         end if
         if (.not. params%enable_species) then
            call fatal_error('input', 'enable_species_enthalpy_diffusion requires transported species')
         end if
         if (params%nspecies <= 0 .and. .not. params%enable_reactions .and. &
             .not. (params%enable_species .and. params%enable_cantera_species)) then
            call fatal_error('input', 'enable_species_enthalpy_diffusion requires nspecies > 0 unless Cantera auto-discovers mechanism species')
         end if
         if (.not. params%enable_cantera_thermo) then
            call fatal_error('input', 'enable_species_enthalpy_diffusion currently requires enable_cantera_thermo=.true. for species sensible enthalpies')
         end if
      end if

      if (params%n_patches < 0 .or. params%n_patches > max_patches) then
         call fatal_error('input', 'n_patches is outside supported range')
      end if

      if (len_trim(params%mesh_dir) == 0) then
         call fatal_error('input', 'mesh_dir cannot be empty')
      end if

      if (len_trim(params%output_dir) == 0) then
         call fatal_error('input', 'output_dir cannot be empty')
      end if

      if (params%vtu_format /= "ascii" .and. params%vtu_format /= "binary") then
         call fatal_error('input', 'vtu_format must be either "ascii" or "binary"')
      end if

      if (params%nspecies < 0 .or. params%nspecies > max_species) then
         call fatal_error('input', 'nspecies is outside supported range')
      end if

      if (params%enable_reactions) then
         if (.not. params%enable_species) call fatal_error('input', 'enable_reactions requires enable_species=.true.')
         if (.not. params%enable_energy) call fatal_error('input', 'enable_reactions requires enable_energy=.true.')
         if (.not. params%enable_cantera_thermo) call fatal_error('input', 'enable_reactions requires enable_cantera_thermo=.true.')
         if (params%chemistry_update_interval <= 0) call fatal_error('input', 'chemistry_update_interval must be positive')
         if (params%chemistry_rtol <= zero) call fatal_error('input', 'chemistry_rtol must be positive')
         if (params%chemistry_atol <= zero) call fatal_error('input', 'chemistry_atol must be positive')
         if (params%chemistry_max_steps <= 0) call fatal_error('input', 'chemistry_max_steps must be positive')
         if (params%chemistry_temperature_cutoff < zero) call fatal_error('input', 'chemistry_temperature_cutoff must be non-negative')
         if (params%chemistry_min_reactive_mass_fraction < zero) call fatal_error('input', 'chemistry_min_reactive_mass_fraction must be non-negative')
         if (params%chemistry_active_species_threshold < zero) call fatal_error('input', 'chemistry_active_species_threshold must be non-negative')
         if (params%chemistry_max_dT_per_step < zero) call fatal_error('input', 'chemistry_max_dT_per_step must be non-negative')
         if (params%chemistry_max_rel_rho_change_per_step < zero) call fatal_error('input', 'chemistry_max_rel_rho_change_per_step must be non-negative')
         if (params%chemistry_max_dY_per_step < zero) call fatal_error('input', 'chemistry_max_dY_per_step must be non-negative')
         if (params%chemistry_dt_safety <= zero .or. params%chemistry_dt_safety > one) call fatal_error('input', 'chemistry_dt_safety must be in (0,1]')
         if (params%chemistry_dt_min_factor <= zero .or. params%chemistry_dt_min_factor > one) call fatal_error('input', 'chemistry_dt_min_factor must be in (0,1]')
         if (params%chemistry_source_relaxation <= zero .or. params%chemistry_source_relaxation > one) call fatal_error('input', 'chemistry_source_relaxation must be in (0,1]')
         if (params%chemistry_n_active_species < 0 .or. params%chemistry_n_active_species > max_species) then
            call fatal_error('input', 'chemistry_n_active_species must be between 0 and max_species')
         end if
      end if

      if (params%variable_density_debug_interval <= 0) then
         call fatal_error('input', 'variable_density_debug_interval must be positive')
      end if

      if (params%enable_ignition_kernel) then
         if (.not. params%enable_energy) call fatal_error('input', 'enable_ignition_kernel requires enable_energy=.true.')
         select case (trim(lowercase(params%ignition_kernel_shape)))
         case ('sphere')
            continue
         case default
            call fatal_error('input', 'ignition_kernel_shape currently supports only sphere')
         end select
         select case (trim(lowercase(params%ignition_kernel_blend)))
         case ('overwrite')
            continue
         case default
            call fatal_error('input', 'ignition_kernel_blend currently supports only overwrite')
         end select
         if (params%ignition_kernel_radius <= zero) call fatal_error('input', 'ignition_kernel_radius must be positive')
         if (params%ignition_kernel_T <= zero) call fatal_error('input', 'ignition_kernel_T must be positive')
         if (len_trim(params%ignition_kernel_composition) > 0 .and. .not. params%enable_species) then
            call fatal_error('input', 'ignition_kernel_composition requires enable_species=.true.')
         end if
      end if

      if (params%restart_from_file .and. len_trim(params%restart_file) == 0) then
         call fatal_error('input', 'restart_from_file=.true. requires restart_file')
      end if
      if (params%write_restart) then
         if (params%restart_interval <= 0) call fatal_error('input', 'write_restart=.true. requires restart_interval > 0')
         if (len_trim(params%restart_output_dir) == 0) call fatal_error('input', 'restart_output_dir cannot be empty')
         if (params%restart_keep < 0) call fatal_error('input', 'restart_keep must be non-negative')
      end if

      if (params%enable_radiation) then
         if (.not. params%enable_energy) call fatal_error('input', 'enable_radiation requires enable_energy=.true.')
         if (params%radiation_update_interval <= 0) call fatal_error('input', 'radiation_update_interval must be positive')
         if (params%radiation_n_wavenumbers <= 0) call fatal_error('input', 'radiation_n_wavenumbers must be positive when radiation is enabled')
         if (params%radiation_n_species < 0 .or. params%radiation_n_species > max_species) then
            call fatal_error('input', 'radiation_n_species is outside supported range')
         end if
         if (params%radiation_n_species > 0 .and. .not. params%enable_species) then
            call fatal_error('input', 'radiation species selection requires enable_species=.true.')
         end if
         select case (trim(lowercase(params%radiation_source_model)))
         case ('none', 'spectral_test', 'external', 'p1', 'dom')
            continue
         case default
            call fatal_error('input', 'radiation_source_model must be one of: none, spectral_test, external, p1, dom')
         end select
         if (trim(lowercase(params%radiation_source_model)) == 'dom') then
            if (trim(lowercase(params%radiation_dom_quadrature)) /= 's2' .and. &
                trim(lowercase(params%radiation_dom_quadrature)) /= 's4' .and. &
                trim(lowercase(params%radiation_dom_quadrature)) /= 's6' .and. &
                trim(lowercase(params%radiation_dom_quadrature)) /= 's8') then
               call fatal_error('input', 'radiation_dom_quadrature must be s2, s4, s6, or s8')
            end if
         end if
         if (params%radiation_n_scalars < 0 .or. params%radiation_n_scalars > max_species) then
            call fatal_error('input', 'radiation_n_scalars is outside supported range')
         end if
         if (trim(lowercase(params%radiation_pressure_source)) /= 'background' .and. &
             trim(lowercase(params%radiation_pressure_source)) /= 'system') then
            call fatal_error('input', 'radiation_pressure_source must be one of: background, system')
         end if
      end if

      call validate_boundary_arrays(params)
   end subroutine validate_params


   subroutine validate_advection_scheme_local(scheme_name, field_name, allow_scalar_high_order)
      character(len=*), intent(in) :: scheme_name
      character(len=*), intent(in) :: field_name
      logical, intent(in) :: allow_scalar_high_order
      character(len=len(scheme_name)) :: scheme

      scheme = trim(lowercase(scheme_name))
      if (len_trim(scheme) == 0) scheme = 'upwind'

      select case (trim(scheme))
      case ('upwind', 'first_order_upwind', 'first-order-upwind', &
            'central', 'central_difference', 'central-difference')
         return
      case ('bounded_central', 'bounded-central', 'bounded_linear', 'bounded-linear', &
            'deferred_central', 'deferred-correction', &
            'muscl', 'limited_linear', 'limited-linear', 'linear_upwind', 'linear-upwind', &
            'higher_order_bounded', 'higher-order-bounded', &
            'quick', 'hlpa', 'quick_limited', 'quick-limited')
         if (allow_scalar_high_order) return
      case default
         continue
      end select

      if (allow_scalar_high_order) then
         call fatal_error('input', trim(field_name)//' must be one of: upwind, central, bounded_central, bounded_linear, muscl, limited_linear, quick, hlpa')
      else
         call fatal_error('input', trim(field_name)//' must be one of: upwind, central')
      end if
   end subroutine validate_advection_scheme_local


   subroutine validate_limiter_name_local(limiter_name)
      character(len=*), intent(in) :: limiter_name
      character(len=len(limiter_name)) :: limiter

      limiter = trim(lowercase(limiter_name))
      if (len_trim(limiter) == 0) limiter = 'barth_jespersen'

      select case (trim(limiter))
      case ('none', 'minmod', 'vanleer', 'van_leer', 'mc', 'monotonized_central', &
            'superbee', 'barth_jespersen', 'barth-jespersen', 'venkatakrishnan', 'hlpa')
         return
      case default
         call fatal_error('input', 'scalar_limiter must be one of: none, minmod, vanleer, mc, superbee, barth_jespersen, venkatakrishnan, hlpa')
      end select
   end subroutine validate_limiter_name_local


   subroutine validate_time_scheme_local(scheme_name, field_name)
      character(len=*), intent(in) :: scheme_name
      character(len=*), intent(in) :: field_name
      character(len=len(scheme_name)) :: scheme

      scheme = trim(lowercase(scheme_name))
      if (len_trim(scheme) == 0) scheme = 'euler'

      select case (trim(scheme))
      case ('euler', 'forward_euler', 'forward-euler', 'ab2', 'adams_bashforth2', 'adams-bashforth2')
         return
      case default
         call fatal_error('input', trim(field_name)//' must be one of: euler, ab2')
      end select
   end subroutine validate_time_scheme_local


   !> Ensures all boundary patch names and types are non-empty.
   subroutine validate_boundary_arrays(params)
      type(case_params_t), intent(in) :: params

      integer :: i

      do i = 1, params%n_patches
         if (len_trim(params%patch_name(i)) == 0) then
            call fatal_error('input', 'patch_name entry cannot be empty')
         end if

         if (len_trim(params%patch_type(i)) == 0) then
            call fatal_error('input', 'patch_type entry cannot be empty')
         end if

         if (params%enable_energy .and. params%patch_T(i) <= zero) then
            call fatal_error('input', 'patch_T entry must be positive when energy is enabled')
         end if
      end do
   end subroutine validate_boundary_arrays

end module mod_input

