mod_mhd_eos.t

Go to the documentation of this file.

!=============================================================================
!> MHD <-> EoS seam: binds the eos% authority into magnetohydrodynamics.
!>
!> mhd_link_eos wires the eos%/phys_ procedure pointers (conversions, pthermal,
!> csound2/gamma1, get_Rfactor, temperature, update_temperature) per eos_type
!> (FI / LTE / PI, with PI energy- and prominence-mode overrides) and per energy
!> formulation (origin / internal_e / hydrodynamic_e / semirelativistic).
!> bind_eos_to_source wires the thermal-conduction / radiative-cooling / thermal-
!> emission / FLD fluid-port callbacks from eos%.
!>
!> Holds the MHD block routines that wrap the shared thermodynamics with MHD's
!> mechanical-energy bookkeeping (kinetic + MAGNETIC energy): the FI/LTE/PI
!> conversions, p_to_e, the pthermal variants, csound2/gamma1, and the PI energy
!> + prominence R-factor / temperature kernels.
!=============================================================================
module mod_mhd_eos
    use mod_global_parameters
    use mod_physics
    use mod_eos
    !> The mode-specific scalar kernels are reached from their sub-modules, not the
    !> mod_eos facade (which no longer re-exports them). Each kernel mpistops if
    !> called under the wrong eos_type/method.
    use mod_eos_lte
    use mod_eos_lte_saha
    use mod_eos_pi
    use mod_eos_container
    use mod_mhd_phys
    use mod_timing
    use mod_radiative_cooling, only: build_y_mod_table
 
    use mod_comm_lib, only: mpistop
 
    implicit none
    private
 
    !> Thermal pressure pointer: set by mhd_link_eos based on energy formulation.
    !> Internal to mod_mhd_eos — external callers use eos%get_thermal_pressure.
    procedure(sub_get_pthermal), pointer, public :: mhd_get_pthermal  => null()
 
    !> Temperature pointer: set by mhd_link_eos based on EoS type and energy formulation.
    procedure(sub_get_pthermal), pointer, public :: mhd_get_temperature => null()
    !> use habitual name of converting to primitive
    procedure(sub_convert), pointer, public :: mhd_to_primitive  => null()
    !> use habitual name of converting to conserved
    procedure(sub_convert), pointer, public :: mhd_to_conserved  => null()
 
    public :: mhd_link_eos
contains
 
    !> Link the appropriate EOS conversion routines based on the selected EoS type
    subroutine mhd_link_eos()
        use mod_usr_methods
 
        !> Primitive <-> Conserved routing
        ! PI (partial approximations) takes the FI conversions as its BASE: in
        ! no-energy mode this is exact (e_int = p/(gamma-1); the variable R-factor
        ! enters only via eos%get_Rfactor, below). In energy mode (eos%ionE) the
        ! PI override block at the end of this routine repoints to_conserved/
        ! to_primitive (and p_to_e/pthermal/csound) at the PI-energy variants, so
        ! the FI assignment here is just a harmless base that gets replaced.
        if (eos%eos_type == 'FI' .or. eos%eos_type == 'PI') then
            ! PI takes its electron count from a variable R factor evaluated on
            ! the rho slot, which under splitting holds only the perturbation.
            if (eos%eos_type == 'PI' .and. has_equi_rho_and_p) &
              call mpistop('PI EoS not supported with equilibrium splitting')
            if(mhd_hydrodynamic_e) then
                eos%to_primitive        => mhd_to_primitive_hde
                eos%to_conserved        => mhd_to_conserved_hde
            else if(mhd_semirelativistic) then
                if(mhd_energy) then
                    eos%to_primitive        => mhd_to_primitive_semirelati
                    eos%to_conserved        => mhd_to_conserved_semirelati
                else
                    eos%to_primitive        => mhd_to_primitive_semirelati_noe
                    eos%to_conserved        => mhd_to_conserved_semirelati_noe
                end if
            else
                if(has_equi_rho_and_p) then
                    eos%to_primitive        => mhd_to_primitive_split_rho
                    eos%to_conserved        => mhd_to_conserved_split_rho
                else if(mhd_internal_e) then
                    eos%to_primitive        => mhd_to_primitive_inte
                    eos%to_conserved        => mhd_to_conserved_inte
                else if(mhd_energy) then
                    eos%to_primitive         => mhd_to_primitive_origin
                    eos%to_conserved         => mhd_to_conserved_origin
                else
                    eos%to_primitive         => mhd_to_primitive_origin_noe
                    eos%to_conserved         => mhd_to_conserved_origin_noe
                end if
            end if
        else if (eos%eos_type == 'LTE') then
            ! Guard unsupported combinations
            if (mhd_hydrodynamic_e) &
              call mpistop('LTE EoS not supported with mhd_hydrodynamic_e')
            if (mhd_semirelativistic) &
              call mpistop('LTE EoS not supported with mhd_semirelativistic')
            if (has_equi_rho_and_p) &
              call mpistop('LTE EoS not supported with equilibrium splitting')
            if (.not. mhd_energy) &
              call mpistop('LTE EoS requires mhd_energy=.true.')
            ! Route LTE conversion based on energy formulation
            if (mhd_internal_e) then
                eos%to_primitive  => mhd_to_primitive_inte_lte
                eos%to_conserved  => mhd_to_conserved_inte_lte
            else
                eos%to_primitive  => mhd_to_primitive_origin_lte
                eos%to_conserved  => mhd_to_conserved_origin_lte
            end if
        else
            call mpistop('Error: Unknown MHD EOS type: ' // trim(eos%eos_type))
        end if
 
        phys_to_primitive       => eos%to_primitive
        phys_to_conserved       => eos%to_conserved
        phys_bind_eos_to_source => bind_eos_to_source
 
        !> p_to_e (pressure -> total energy for origin, or eint for inte)
        if (mhd_internal_e) then
            eos%p_to_e => mhd_p_to_eint
        else
            eos%p_to_e => mhd_p_to_e
        end if
 
        !> Sound speed and Gamma1
        if (eos%eos_type == 'LTE' .and. eos%ionE) then
            eos%get_csound2 => mhd_get_csound2_lte
            phys_get_gamma1 => mhd_get_gamma1_lte
        else
            eos%get_csound2 => mhd_get_csound2_fi
            phys_get_gamma1 => get_gamma1_fi
        end if
 
        !> Prolongation
        if (eos%eos_type == 'LTE' .and. eos%ionE) then
            phys_to_prolong   => mhd_to_prolong_lte
            phys_from_prolong => mhd_from_prolong_lte
        end if
 
        !> Rfactor: only the FI case is physics-dependent (the usr_Rfactor user
        !> hook). The LTE and PI bindings are pure eos% routines, so they are set
        !> in eos_finalise_{LTE,PI} (which run after this and before
        !> bind_eos_to_source, so they win) -- keeping those targets private.
        if(associated(usr_rfactor)) then
            eos%get_Rfactor => usr_rfactor
        else
            eos%get_Rfactor => rfactor_from_constant_ionization
        end if
 
 
        !> Thermal pressure routing
        if(mhd_internal_e) then
            phys_get_pthermal        => mhd_get_pthermal_inte
            mhd_get_pthermal         => mhd_get_pthermal_inte
        else if(mhd_hydrodynamic_e) then
            phys_get_pthermal        => mhd_get_pthermal_hde
            mhd_get_pthermal         => mhd_get_pthermal_hde
        else if(mhd_semirelativistic) then
            phys_get_pthermal        => mhd_get_pthermal_semirelati
            mhd_get_pthermal         => mhd_get_pthermal_semirelati
        else if(mhd_energy) then
            phys_get_pthermal        => mhd_get_pthermal_origin
            mhd_get_pthermal         => mhd_get_pthermal_origin
        else
            phys_get_pthermal        => mhd_get_pthermal_noe
            mhd_get_pthermal         => mhd_get_pthermal_noe
        end if
 
        ! For LTE, override with version that delegates to eos%get_thermal_pressure
        if (eos%eos_type == 'LTE') then
            phys_get_pthermal => mhd_get_pthermal_lte
            mhd_get_pthermal  => mhd_get_pthermal_lte
        end if
        ! Single entry point: all callers use eos%get_thermal_pressure
        eos%get_thermal_pressure => mhd_get_pthermal
 
        !> Temperature routing
        if(eos%eos_type == 'PI' .or. eos%eos_type == 'LTE') then
            mhd_get_temperature => mhd_get_temperature_from_te
        else
            if(mhd_internal_e) then
                if(has_equi_rho_and_p) then
                    mhd_get_temperature => mhd_get_temperature_from_eint_with_equi
                else
                    mhd_get_temperature => mhd_get_temperature_from_eint
                end if
            else
                if(has_equi_rho_and_p) then
                    mhd_get_temperature => mhd_get_temperature_from_etot_with_equi
                else
                    mhd_get_temperature => mhd_get_temperature_from_etot
                end if
            end if
        end if
 
        !> Internal energy extraction
        if(mhd_internal_e) then
            phys_get_ei => mhd_get_ei_inte
        else if(mhd_energy) then
            phys_get_ei => mhd_get_ei_origin
        end if
 
        !> PI energy-mode overrides (eos_type='PI', ionE=.true.)
        ! PI shares FI's normalisation and no-energy path verbatim; energy mode
        ! differs ONLY in the eint<->p relation (eint carries the ionisation
        ! energy). Override exactly the routines that touch that relation, all
        ! delegating to the portable scalar backend in mod_eos_PI. FI
        ! routines stay pristine (no flag branch in the hot path).
        if (eos%eos_type == 'PI' .and. eos%ionE) then
            if (.not. mhd_energy) &
              call mpistop('PI energy EoS requires mhd_energy=.true.')
            if (mhd_hydrodynamic_e) &
              call mpistop('PI energy EoS not supported with mhd_hydrodynamic_e')
            if (mhd_semirelativistic) &
              call mpistop('PI energy EoS not supported with mhd_semirelativistic')
            if (has_equi_rho_and_p) &
              call mpistop('PI energy EoS not supported with equilibrium splitting')
 
            if (mhd_internal_e) then
                eos%to_conserved  => mhd_to_conserved_inte_pi
                eos%to_primitive  => mhd_to_primitive_inte_pi
                eos%p_to_e        => mhd_p_to_eint_pi
                phys_get_pthermal => mhd_get_pthermal_inte_pi
                mhd_get_pthermal  => mhd_get_pthermal_inte_pi
            else
                eos%to_conserved  => mhd_to_conserved_origin_pi
                eos%to_primitive  => mhd_to_primitive_origin_pi
                eos%p_to_e        => mhd_p_to_e_pi
                phys_get_pthermal => mhd_get_pthermal_origin_pi
                mhd_get_pthermal  => mhd_get_pthermal_origin_pi
            end if
            phys_to_primitive        => eos%to_primitive
            phys_to_conserved        => eos%to_conserved
            eos%get_thermal_pressure => mhd_get_pthermal
 
            !> eos%get_csound2 => get_csound2_PI is physics-independent: set in
            !> eos_finalise_PI (ionE arm) so the target stays private.
            phys_get_gamma1 => get_gamma1_pi
 
            ! Te_ refreshed from eint each substep (direct inversion, no lag).
            ! Module-specific (uses MHD's Te_ index): PI registers Te via
            ! var_set_auxvar, which does NOT set the generic iw_te, so the
            ! temperature aux must be addressed through the module's own Te_.
        end if
 
        ! phys_e_to_ei / phys_ei_to_e assigned in mhd_phys_init (mod_mhd_phys.t)
        mhd_to_primitive => eos%to_primitive
        mhd_to_conserved => eos%to_conserved
 
    end subroutine mhd_link_eos
 
    !> Called by eos_finalise via phys_bind_eos_to_source to link
    !> TC and RC modules to EoS-aware function pointers.
    subroutine bind_eos_to_source()
 
        ! Override eos%get_temperature_from_{etot,eint} per (eos_type, internal_e, equi).
        ! Runs AFTER eos_finalise's unconditional generic-helper assignment.
        ! Two semantically distinct hooks:
        !   - from_etot: caller's w(:,e_) is total energy; subtract KE+ME before T.
        !   - from_eint: caller's w(:,e_) is gas internal energy; compute T directly.
        ! In mhd_internal_e mode w(:,e_) IS already eint, so the etot hook is bound
        ! to the same routine as the eint hook (no subtraction). Required because
        ! eos_finalise's generic helper calls phys_e_to_ei, which is left unbound
        ! by mhd_phys_init when mhd_internal_e=.true.
        if (eos%eos_type == 'LTE') then
            if (mhd_internal_e) then
                ! w(:,e_) is eint; both hooks do the LTE table lookup directly.
                eos%get_temperature_from_etot => eos%get_temperature_from_eint
            else
                ! Total-energy: subtract KE+ME (via phys_e_to_ei) then LTE table lookup.
                eos%get_temperature_from_etot => mhd_get_temperature_from_etot_lte
            end if
        else  ! FI
            if (mhd_internal_e) then
                ! w(:,e_) is eint; both hooks share the eint variant.
                if (has_equi_rho_and_p) then
                    eos%get_temperature_from_etot => mhd_get_temperature_from_eint_with_equi
                    eos%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
                else
                    eos%get_temperature_from_etot => mhd_get_temperature_from_eint
                    eos%get_temperature_from_eint => mhd_get_temperature_from_eint
                end if
            else
                ! Total-energy: distinct etot/eint paths.
                if (has_equi_rho_and_p) then
                    eos%get_temperature_from_etot => mhd_get_temperature_from_etot_with_equi
                    eos%get_temperature_from_eint => mhd_get_temperature_from_eint_with_equi
                else
                    eos%get_temperature_from_etot => mhd_get_temperature_from_etot
                    eos%get_temperature_from_eint => mhd_get_temperature_from_eint
                end if
            end if
        end if
 
        if (allocated(tc_fl)) then
            tc_fl%get_temperature_from_conserved => eos%get_temperature_from_etot
            if (eos%eos_type == 'LTE' .and. eos%ionE) then
                tc_fl%get_temperature_from_eint => get_temperature_from_eint_fast_lte
            else
                tc_fl%get_temperature_from_eint => eos%get_temperature_from_eint
            end if
            tc_fl%get_rho => eos%get_rho
            tc_fl%get_ne_nH => eos%get_ne_nH
            tc_fl%get_var_Rfactor => eos%get_Rfactor
            tc_fl%inv_gamma_minus_1 =  eos%inv_gamma_minus_1
            tc_fl%nH2rhoFactor      =  eos%nH2rhoFactor
            tc_fl%log_T_floor       =  eos_get_log_t_floor()
            tc_fl%eint_from_T       => eint_nh_from_t
            ! Equilibrium-specific pointers
            if(has_equi_rho_and_p .and. mhd_equi_thermal) then
                tc_fl%subtract_equi = .true.
                tc_fl%get_temperature_equi => mhd_get_temperature_equi
                tc_fl%get_rho_equi => mhd_get_rho_equi
            else
                tc_fl%subtract_equi = .false.
            end if
        end if
 
        if (allocated(rc_fl)) then
            rc_fl%get_rho          => eos%get_rho
            rc_fl%get_pthermal     => eos%get_thermal_pressure
            rc_fl%get_var_Rfactor  => eos%get_Rfactor
            rc_fl%get_Te           => eos%get_Te
            rc_fl%get_ne_nH        => eos%get_ne_nH
            nullify(rc_fl%get_rho2_factor)
            if(mhd_uawsom) rc_fl%get_rho2_factor => mhd_uawsom_rho2_factor
            rc_fl%ionE              =  eos%ionE
            rc_fl%method            =  eos%method
            rc_fl%inv_gamma_minus_1 =  eos%inv_gamma_minus_1
            rc_fl%nH2rhoFactor      =  eos%nH2rhoFactor
            rc_fl%eion_per_nH       =  eos%eion_per_nH
            rc_fl%eint_from_T       => eint_nh_from_t
            rc_fl%p2eint            => p2eint_from_nh_p
            rc_fl%T_from_eint       => t_from_nh_eint
            rc_fl%y_from_eint       => y_from_nh_eint
            ! Equilibrium-specific pointers
            if(has_equi_rho_and_p .and. mhd_equi_thermal) then
                rc_fl%subtract_equi = .true.
                rc_fl%get_rho_equi => mhd_get_rho_equi
                rc_fl%get_pthermal_equi => mhd_get_pe_equi
                rc_fl%get_ne_nH_equi => mhd_get_ne_nh_equi
                rc_fl%get_temperature_equi => mhd_get_temperature_equi
            else
                rc_fl%subtract_equi = .false.
            end if
            !> Build the variable-c_V Townsend Y_mod table now that all
            !> EoS tables (eint_from_T, T, neOnH) are in code units.
            !> build_Y_mod_table checks coolmethod=='exact' and .not.isPPL
            !> internally and early-returns otherwise. Feed it the inverse-table nH
            !> grid via the port so it never touches eos% directly. LTE only:
            !> eos_get_eintT_grid reads the LTE eintT grid, and PI (T-only) cooling
            !> runs the classical Townsend path (Y_mod-for-PI is a later refinement).
            if (eos%ionE .and. eos%eos_type == 'LTE') then
                call eos_get_eintt_grid(rc_fl%Y_mod_n_nH, &
                     rc_fl%Y_mod_lg_nH_min, rc_fl%Y_mod_lg_nH_max)
                call build_y_mod_table(rc_fl)
            end if
        end if
 
        !> PI energy mode: the eint<->T scalar callbacks wired above are the
        !> LTE-table lookups, which PI does not load. Repoint cooling/conduction
        !> at the ionisation backend (rho=nH=1 per H). Classical Townsend; Y_mod
        !> is not built for PI (guard above), so the classical branch is used.
        if (eos%eos_type == 'PI' .and. eos%ionE) then
            if (allocated(tc_fl)) tc_fl%eint_from_T => eint_from_t_pi
            if (allocated(rc_fl)) then
                rc_fl%eint_from_T => eint_from_t_pi
                rc_fl%p2eint      => p2eint_pi
                rc_fl%T_from_eint => t_from_eint_pi
                rc_fl%y_from_eint => y_from_eint_pi
            end if
        end if
 
        if (allocated(te_fl_mhd)) then
            te_fl_mhd%get_rho          => eos%get_rho
            te_fl_mhd%get_pthermal     => eos%get_thermal_pressure
            te_fl_mhd%get_var_Rfactor  => eos%get_Rfactor
            te_fl_mhd%get_ne_nH        => eos%get_ne_nH
        end if
 
        if (allocated(fld_fl)) then
            !> Radiation (FLD) fluid: gas-EoS callbacks. get_temperature_from_pressure
            !> is the FI T=p/(R*rho) routine (the LTE variant is a later pass).
            fld_fl%gamma       =  eos%gamma
            fld_fl%get_tgas    => eos%get_temperature_from_pressure
            fld_fl%get_Rfactor => eos%get_Rfactor
        end if
 
    end subroutine bind_eos_to_source
 
    !> FI conversion routines (signatures updated for convert_condition interface)
 
    !> Transform primitive variables into conservative ones (origin, total energy)
    subroutine mhd_to_conserved_origin(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! Calculate total energy from pressure, kinetic and magnetic energy
            w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1&
                      +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                      +(^c&w(ix^d,b^c_)**2+))
            if(mhd_uawsom) w(ix^d,e_)=w(ix^d,e_)+&
                 w(ix^d,waplus_)+w(ix^d,waminus_)+w(ix^d,wkplus_)+w(ix^d,wkminus_)
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_to_conserved_origin
 
    !> Transform primitive variables into conservative ones (no energy)
    subroutine mhd_to_conserved_origin_noe(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_to_conserved_origin_noe
 
    !> Transform primitive variables into conservative ones (hydrodynamic energy)
    subroutine mhd_to_conserved_hde(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! Calculate total energy from pressure, kinetic and magnetic energy
            w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1&
                      +half*(^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_to_conserved_hde
 
    !> Transform primitive variables into conservative ones (internal energy)
    subroutine mhd_to_conserved_inte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! Calculate internal energy from pressure
            w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_to_conserved_inte
 
    !> Transform primitive variables into conservative ones (split rho)
    subroutine mhd_to_conserved_split_rho(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: rho
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            rho=w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i)
            ! Calculate total energy from pressure, kinetic and magnetic energy
            w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1&
                      +half*((^c&w(ix^d,m^c_)**2+)*rho&
                            +(^c&w(ix^d,b^c_)**2+))
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=rho*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_to_conserved_split_rho
 
    !> Transform primitive variables into conservative ones (semirelativistic)
    subroutine mhd_to_conserved_semirelati(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: e(ixo^s,1:ndir), s(ixo^s,1:ndir)
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            {^ifthreec
                e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
                e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
                e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                s(ix^d,1)=e(ix^d,2)*w(ix^d,b3_)-e(ix^d,3)*w(ix^d,b2_)
                s(ix^d,2)=e(ix^d,3)*w(ix^d,b1_)-e(ix^d,1)*w(ix^d,b3_)
                s(ix^d,3)=e(ix^d,1)*w(ix^d,b2_)-e(ix^d,2)*w(ix^d,b1_)
            }
            {^iftwoc
                e(ix^d,1)=zero
                ! switch 3 with 2 to add 3 when ^C from 1 to 2
                e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                s(ix^d,1)=-e(ix^d,2)*w(ix^d,b2_)
                s(ix^d,2)=e(ix^d,2)*w(ix^d,b1_)
            }
            {^ifonec
                e(ix^d,1)=zero
                s(ix^d,1)=zero
            }
            if(mhd_internal_e) then
                ! internal energy
                w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1
            else
                ! equation (9)
                ! Calculate total energy from internal, kinetic and magnetic energy
                w(ix^d,e_)=w(ix^d,p_)*eos%inv_gamma_minus_1&
                          +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                          +(^c&w(ix^d,b^c_)**2+)&
                          +(^c&e(ix^d,^c)**2+)*eos%inv_squared_c)
            end if
 
            ! Convert velocity to momentum, equation (9)
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)+s(ix^d,^c)*eos%inv_squared_c\
 
        {end do\}
 
    end subroutine mhd_to_conserved_semirelati
 
    subroutine mhd_to_conserved_semirelati_noe(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: e(ixo^s,1:ndir), s(ixo^s,1:ndir)
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            {^ifthreec
                e(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
                e(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
                e(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                s(ix^d,1)=e(ix^d,2)*w(ix^d,b3_)-e(ix^d,3)*w(ix^d,b2_)
                s(ix^d,2)=e(ix^d,3)*w(ix^d,b1_)-e(ix^d,1)*w(ix^d,b3_)
                s(ix^d,3)=e(ix^d,1)*w(ix^d,b2_)-e(ix^d,2)*w(ix^d,b1_)
            }
            {^iftwoc
                e(ix^d,1)=zero
                ! switch 3 with 2 to add 3 when ^C from 1 to 2
                e(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                s(ix^d,1)=-e(ix^d,2)*w(ix^d,b2_)
                s(ix^d,2)=e(ix^d,2)*w(ix^d,b1_)
            }
            {^ifonec
                s(ix^d,1)=zero
            }
            ! Convert velocity to momentum, equation (9)
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)+s(ix^d,^c)*eos%inv_squared_c\
 
        {end do\}
 
    end subroutine mhd_to_conserved_semirelati_noe
 
    !> Transform conservative variables into primitive ones (origin)
    subroutine mhd_to_primitive_origin(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_origin')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Calculate pressure = (gamma-1) * (e-ek-eb)
            w(ix^d,p_)=eos%gamma_minus_1*(w(ix^d,e_)&
                      -half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+)&
                        +(^c&w(ix^d,b^c_)**2+))&
                      -mhd_uawsom_wave_energy_cell(w(ix^d,:)))
        {end do\}
 
    end subroutine mhd_to_primitive_origin
 
    !> Transform conservative variables into primitive ones (no energy)
    subroutine mhd_to_primitive_origin_noe(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_origin_noe')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
        {end do\}
 
    end subroutine mhd_to_primitive_origin_noe
 
    !> Transform conservative variables into primitive ones (hde)
    subroutine mhd_to_primitive_hde(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        integer                         :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_hde')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Calculate pressure = (gamma-1) * (e-ek)
            w(ix^d,p_)=eos%gamma_minus_1*(w(ix^d,e_)-half*w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+))
        {end do\}
 
    end subroutine mhd_to_primitive_hde
 
    !> Transform conservative variables into primitive ones (internal energy)
    subroutine mhd_to_primitive_inte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision                :: inv_rho
        integer                         :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_inte')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! Calculate pressure = (gamma-1) * e_internal
            w(ix^d,p_)=w(ix^d,e_)*eos%gamma_minus_1
            ! Convert momentum to velocity
            inv_rho = 1.d0/w(ix^d,rho_)
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
        {end do\}
 
    end subroutine mhd_to_primitive_inte
 
    !> Transform conservative variables into primitive ones (split rho)
    subroutine mhd_to_primitive_split_rho(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_split_rho')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho=1.d0/(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Calculate pressure = (gamma-1) * (e-ek-eb)
            w(ix^d,p_)=eos%gamma_minus_1*(w(ix^d,e_)&
                        -half*((w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,b0i))*&
                        (^c&w(ix^d,m^c_)**2+)+(^c&w(ix^d,b^c_)**2+)))
        {end do\}
 
    end subroutine mhd_to_primitive_split_rho
 
    !> Transform conservative variables into primitive ones (semirelativistic)
    subroutine mhd_to_primitive_semirelati(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: b(ixo^s,1:ndir), tmp, b2, gamma2, inv_rho
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_semirelati')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            b2=(^c&w(ix^d,b^c_)**2+)
            if(b2>smalldouble) then
                tmp=1.d0/sqrt(b2)
            else
                tmp=0.d0
            end if
            ^c&b(ix^d,^c)=w(ix^d,b^c_)*tmp\
            tmp=(^c&b(ix^d,^c)*w(ix^d,m^c_)+)
 
            inv_rho=1.d0/w(ix^d,rho_)
            ! Va^2/c^2
            b2=b2*inv_rho*eos%inv_squared_c
            ! equation (15)
            gamma2=1.d0/(1.d0+b2)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=gamma2*(w(ix^d,m^c_)+b2*b(ix^d,^c)*tmp)*inv_rho\
 
            if(mhd_internal_e) then
                ! internal energy to pressure
                w(ix^d,p_)=eos%gamma_minus_1*w(ix^d,e_)
            else
                ! E=Bxv
                {^ifthreec
                    b(ix^d,1)=w(ix^d,b2_)*w(ix^d,m3_)-w(ix^d,b3_)*w(ix^d,m2_)
                    b(ix^d,2)=w(ix^d,b3_)*w(ix^d,m1_)-w(ix^d,b1_)*w(ix^d,m3_)
                    b(ix^d,3)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                }
                {^iftwoc
                    b(ix^d,1)=zero
                    b(ix^d,2)=w(ix^d,b1_)*w(ix^d,m2_)-w(ix^d,b2_)*w(ix^d,m1_)
                }
                {^ifonec
                    b(ix^d,1)=zero
                }
                ! Calculate pressure = (gamma-1) * (e-eK-eB-eE)
                w(ix^d,p_)=eos%gamma_minus_1*(w(ix^d,e_)&
                          -half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                          +(^c&w(ix^d,b^c_)**2+)&
                          +(^c&b(ix^d,^c)**2+)*eos%inv_squared_c))
            end if
        {end do\}
 
    end subroutine mhd_to_primitive_semirelati
 
    !> Transform conservative variables into primitive ones (semirelativistic noe)
    subroutine mhd_to_primitive_semirelati_noe(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: b(ixo^s,1:ndir),tmp,b2,gamma2,inv_rho
        integer :: ix^d, idir
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, 'mhd_to_primitive_semirelati_noe')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            b2=(^c&w(ix^d,b^c_)**2+)
            if(b2>smalldouble) then
                tmp=1.d0/sqrt(b2)
            else
                tmp=0.d0
            end if
            ^c&b(ix^d,^c)=w(ix^d,b^c_)*tmp\
            tmp=(^c&b(ix^d,^c)*w(ix^d,m^c_)+)
 
            inv_rho=1.d0/w(ix^d,rho_)
            ! Va^2/c^2
            b2=b2*inv_rho*eos%inv_squared_c
            ! equation (15)
            gamma2=1.d0/(1.d0+b2)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=gamma2*(w(ix^d,m^c_)+b2*b(ix^d,^c)*tmp)*inv_rho\
        {end do\}
 
    end subroutine mhd_to_primitive_semirelati_noe
 
    !> LTE conversion routines
 
    !> LTE: Primitive (rho, v, p, B) -> Conserved (rho, rho*v, E_total, B)
    !> E_total = eint + 0.5*rho*v^2 + 0.5*B^2
    subroutine mhd_to_conserved_origin_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        timeeos0 = mpi_wtime()
 
        ! Convert p -> E_total (eint + eK + eB) via EoS table
        call mhd_p_to_e(ixi^l, ixo^l, w, x)
        ! Convert velocity to momentum
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine mhd_to_conserved_origin_lte
 
    !> LTE: Primitive (rho, v, p, B) -> Conserved (rho, rho*v, eint, B)
    !> Internal energy formulation: e_ stores eint only
    subroutine mhd_to_conserved_inte_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        timeeos0 = mpi_wtime()
 
        ! Convert p -> eint via EoS table (no mechanical energy)
        call mhd_p_to_eint(ixi^l, ixo^l, w, x)
        ! Convert velocity to momentum
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine mhd_to_conserved_inte_lte
 
    !> Convert pressure to total energy (eint + eK + eB) for origin formulation.
    !> Uses p2eint table for LTE ionE, with FI bypass for hot cells.
    subroutine mhd_p_to_e(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
        double precision :: p_to_eint, p_over_rho
        double precision :: nh(ixi^s), nh_in(ixi^s), p_in(ixi^s)
        double precision :: t_solve, y_solve, eint_nh_solve
        double precision :: log_eint_mid, eint_total
 
        if (eos%ionE) then
            call eos%get_nH(w, x, ixi^l, ixo^l, nh)
            nh_in(ixo^s) = dlog10(nh(ixo^s))
            if (eos%p2eint_method /= 'bisect') then
                p_in(ixo^s) = dlog10(w(ixo^s,p_)) - nh_in(ixo^s)
            end if
        endif
 
        p_to_eint = eos%inv_gamma_minus_1
        {do ix^db=ixomin^db,ixomax^db\}
            if (eos%ionE) then
                p_over_rho = w(ix^d,p_) / w(ix^d,rho_)
                if (p_over_rho > eos%p_rho_FI_threshold) then
                    p_to_eint = eos%inv_gamma_minus_1 &
                        + eos%eion_per_nH * nh(ix^d) / w(ix^d,p_)
                    w(ix^d,e_)=w(ix^d,p_)*p_to_eint&
                              +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                              +(^c&w(ix^d,b^c_)**2+))
                else if (eos%method == 'analytic') then
                    call saha_state_from_nh_p(nh(ix^d), w(ix^d,p_), &
                        t_solve, y_solve, eint_nh_solve)
                    w(ix^d,e_) = eint_nh_solve * nh(ix^d) &
                        +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                        +(^c&w(ix^d,b^c_)**2+))
                else if (eos%p2eint_method == 'bisect') then
                    call eint_from_p_bisect(nh_in(ix^d), &
                        dlog10(w(ix^d,p_)), log_eint_mid)
                    eint_total = nh(ix^d) * 10.0d0**log_eint_mid
                    eint_total = max(eint_total, &
                        nh(ix^d) * 10.0d0**eos%T%var2_min)
                    w(ix^d,e_) = eint_total + &
                        half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                        +(^c&w(ix^d,b^c_)**2+))
                else
                    p_to_eint = p2eint_from_nh_p(nh_in(ix^d), p_in(ix^d))
                    w(ix^d,e_)=w(ix^d,p_)*p_to_eint&
                              +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                              +(^c&w(ix^d,b^c_)**2+))
                end if
            else
                w(ix^d,e_)=w(ix^d,p_)*p_to_eint&
                          +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                          +(^c&w(ix^d,b^c_)**2+))
            end if
            if(mhd_uawsom) w(ix^d,e_)=w(ix^d,e_)+&
                 w(ix^d,waplus_)+w(ix^d,waminus_)+w(ix^d,wkplus_)+w(ix^d,wkminus_)
        {end do\}
 
    end subroutine mhd_p_to_e
 
    !> Convert pressure to internal energy only (for internal_e formulation).
    !> Uses p2eint table for LTE ionE, with FI bypass for hot cells.
    subroutine mhd_p_to_eint(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        integer :: ix^d
        double precision :: p_to_eint, p_over_rho
        double precision :: nh(ixi^s), nh_in(ixi^s), p_in(ixi^s)
        double precision :: t_solve, y_solve, eint_nh_solve
        double precision :: log_eint_mid, eint_total
 
        if (eos%ionE) then
            call eos%get_nH(w, x, ixi^l, ixo^l, nh)
            nh_in(ixo^s) = dlog10(nh(ixo^s))
            if (eos%p2eint_method /= 'bisect') then
                p_in(ixo^s) = dlog10(w(ixo^s,p_)) - nh_in(ixo^s)
            end if
        endif
 
        p_to_eint = eos%inv_gamma_minus_1
        {do ix^db=ixomin^db,ixomax^db\}
            if (eos%ionE) then
                p_over_rho = w(ix^d,p_) / w(ix^d,rho_)
                if (p_over_rho > eos%p_rho_FI_threshold) then
                    p_to_eint = eos%inv_gamma_minus_1 &
                        + eos%eion_per_nH * nh(ix^d) / w(ix^d,p_)
                    w(ix^d,e_) = w(ix^d,p_) * p_to_eint
                else if (eos%method == 'analytic') then
                    call saha_state_from_nh_p(nh(ix^d), w(ix^d,p_), &
                        t_solve, y_solve, eint_nh_solve)
                    w(ix^d,e_) = eint_nh_solve * nh(ix^d)
                else if (eos%p2eint_method == 'bisect') then
                    call eint_from_p_bisect(nh_in(ix^d), &
                        dlog10(w(ix^d,p_)), log_eint_mid)
                    eint_total = nh(ix^d) * 10.0d0**log_eint_mid
                    eint_total = max(eint_total, &
                        nh(ix^d) * 10.0d0**eos%T%var2_min)
                    w(ix^d,e_) = eint_total
                else
                    p_to_eint = p2eint_from_nh_p(nh_in(ix^d), p_in(ix^d))
                    w(ix^d,e_) = w(ix^d,p_) * p_to_eint
                end if
            else
                w(ix^d,e_) = w(ix^d,p_) * p_to_eint
            end if
        {end do\}
 
    end subroutine mhd_p_to_eint
 
    !> LTE: Conserved (rho, rho*v, E_total, B) -> Primitive (rho, v, p, B)
    !> Uses T and neOnH tables with FI bypass for hot cells.
    subroutine mhd_to_primitive_origin_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, eint_in
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'mhd_to_primitive_origin_LTE')
        end if
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        if (eos%ionE) then
            log_nh(ixo^s) = dlog10(nh(ixo^s))
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! Extract internal energy: eint = E - 0.5*rho*v^2 - 0.5*B^2
            eint_val = w(ix^d,e_) &
                      - half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+) &
                      + (^c&w(ix^d,b^c_)**2+))
            ! Floor eint to prevent unphysical values
            if (eos%method /= 'analytic') then
                eint_val = max(eint_val, nh(ix^d) * 10.0d0**eos%T%var2_min)
            end if
            eint_val = max(eint_val, smalldouble)
 
            if (eos%ionE) then
                if (eint_val * inv_rho > eos%eint_rho_FI_threshold) then
                    ! FI bypass: p = (gamma-1) * (eint - eion*nH)
                    w(ix^d,p_) = eos%gamma_minus_1 &
                        * (eint_val - eos%eion_per_nH * nh(ix^d))
                else
                    ! Ionisation zone: single p/nH lookup (replaces separate T + y)
                    eint_in = dlog10(eint_val) - log_nh(ix^d)
                    w(ix^d,p_) = nh(ix^d) * p_nh_from_eint(log_nh(ix^d), eint_in)
                end if
            else
                ! FI without ionE tables
                w(ix^d,p_) = eos%gamma_minus_1 * eint_val
            end if
        {end do\}
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine mhd_to_primitive_origin_lte
 
    !> LTE: Conserved (rho, rho*v, eint, B) -> Primitive (rho, v, p, B)
    !> Internal energy formulation: eint = w(e_) directly.
    subroutine mhd_to_primitive_inte_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, eint_in, t_loc, y_loc
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'mhd_to_primitive_inte_LTE')
        end if
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        if (eos%ionE) then
            log_nh(ixo^s) = dlog10(nh(ixo^s))
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! eint = w(e_) directly for internal energy formulation
            eint_val = w(ix^d,e_)
            eint_val = max(eint_val, nh(ix^d) * 10.0d0**eos%T%var2_min)
 
            if (eos%ionE) then
                if (eint_val / w(ix^d,rho_) > eos%eint_rho_FI_threshold) then
                    w(ix^d,p_) = eos%gamma_minus_1 &
                        * (eint_val - eos%eion_per_nH * nh(ix^d))
                else
                    eint_in = dlog10(eint_val) - log_nh(ix^d)
                    w(ix^d,p_) = nh(ix^d) * p_nh_from_eint(log_nh(ix^d), eint_in)
                end if
            else
                w(ix^d,p_) = eos%gamma_minus_1 * eint_val
            end if
 
            ! Convert momentum to velocity
            inv_rho = 1.d0/w(ix^d,rho_)
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
        {end do\}
 
        timeeos_conv=timeeos_conv+(mpi_wtime()-timeeos0)
 
    end subroutine mhd_to_primitive_inte_lte
 
    !> LTE prolongation routines
 
    !> Convert conserved (rho, rho*v, E_total, B) to prolongation form (rho, v, T, B).
    !> Temperature is stored in the p_ slot for AMR interpolation.
    !> T is smoothest through the ionisation zone (conduction-dominated, linear gradient).
    subroutine mhd_to_prolong_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, t_loc, y_loc
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            if (mhd_internal_e) then
                eint_val = w(ix^d,e_)
            else
                ! Compute eint = E - 0.5*rho*v^2 - 0.5*B^2
                eint_val = w(ix^d,e_) &
                    - half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+) &
                    + (^c&w(ix^d,b^c_)**2+))
            end if
            ! Floor eint to prevent unphysical values
            if (eos%method /= 'analytic') then
                eint_val = max(eint_val, nh(ix^d) * 10.0d0**eos%T%var2_min)
            end if
            eint_val = max(eint_val, smalldouble)
 
            if (eint_val * inv_rho > eos%eint_rho_FI_threshold) then
                ! FI: T = (gamma-1)*(eint - eion*nH) / (nH * n_per_nH_FI)
                w(ix^d,p_) = eos%gamma_minus_1 &
                    * (eint_val - eos%eion_per_nH * nh(ix^d)) &
                    / (nh(ix^d) * eos%n_per_nH_FI)
            else if (eos%method == 'analytic') then
                ! Analytical Saha: solve for T from eint
                call saha_t_from_nh_eint(nh(ix^d), &
                    eint_val / nh(ix^d), t_loc, y_loc)
                w(ix^d,p_) = t_loc
            else
                ! Ionisation zone: T from table
                w(ix^d,p_) = t_from_nh_eint( &
                    log_nh(ix^d), &
                    dlog10(eint_val) - log_nh(ix^d))
            end if
        {end do\}
 
    end subroutine mhd_to_prolong_lte
 
    !> Convert prolongation form (rho, v, T, B) to conserved (rho, rho*v, E_total, B).
    !> T is read from the p_ slot. Uses eint_nH_from_T table for back-conversion.
    subroutine mhd_from_prolong_lte(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: t_val, eint_val, t_fi, log_t_min
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        ! Temperature above which gas is fully ionised
        t_fi = (eos%eint_rho_FI_threshold &
            * eos%nH2rhoFactor - eos%eion_per_nH) &
            * eos%gamma_minus_1 / eos%n_per_nH_FI
 
        ! Floor for log_T into the (rho, T) inverse table; entropy method uses
        ! eos%eintT, legacy 'tables' method uses eos%eint_from_T. Picking the
        ! wrong container leaves var2_min = 0 → floors T at 10^6 K.
        if (eos%method == 'entropy') then
            log_t_min = eos%eintT%var2_min
        else
            log_t_min = eos%eint_from_T%var2_min
        end if
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            t_val = w(ix^d,p_)  ! T stored in p_ slot
            if (t_val > t_fi) then
                ! FI: eint = nH*(n_per_nH*T/(gamma-1) + eion)
                eint_val = nh(ix^d) &
                    * (eos%n_per_nH_FI * t_val * eos%inv_gamma_minus_1 &
                    + eos%eion_per_nH)
            else if (eos%method == 'analytic') then
                ! Analytical Saha: eint from T directly
                eint_val = saha_eint_from_nh_t(nh(ix^d), t_val) * nh(ix^d)
            else
                ! Ionisation zone: eint/nH from T table
                eint_val = eint_nh_from_t( &
                    log_nh(ix^d), &
                    dlog10(max(t_val, 10.0d0**log_t_min))) &
                    * nh(ix^d)
            end if
            if (mhd_internal_e) then
                w(ix^d,e_) = eint_val
            else
                ! E = eint + 0.5*rho*v^2 + 0.5*B^2
                w(ix^d,e_) = eint_val &
                    + half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+) &
                    + (^c&w(ix^d,b^c_)**2+))
            end if
            ! Convert velocity to momentum
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)\
        {end do\}
 
    end subroutine mhd_from_prolong_lte
 
    !> Sound speed routines
 
    !> Acoustic sound speed squared for FI (constant gamma) EoS.
    !> Expects w in primitive form: w(p_) = pressure, w(rho_) = density.
    subroutine mhd_get_csound2_fi(w, x, ixI^L, ixO^L, cs2)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: cs2(ixi^s)
 
        double precision :: rho(ixi^s), pth(ixi^s)
 
        timeeos0 = mpi_wtime()
 
        !> Both slots hold perturbations under equilibrium splitting, so the
        !> background is restored before the ratio is formed.
        call eos%get_rho(w, x, ixi^l, ixo^l, rho)
        pth(ixo^s) = w(ixo^s, p_)
        if (iw_equi_p > 0) pth(ixo^s) = pth(ixo^s) &
             + block%equi_vars(ixo^s, iw_equi_p, b0i)
        cs2(ixo^s) = eos%gamma * pth(ixo^s) / rho(ixo^s)
 
        timeeos_csound = timeeos_csound + (mpi_wtime()-timeeos0)
 
    end subroutine mhd_get_csound2_fi
 
    !> Acoustic sound speed squared for LTE+IonE EoS using Gamma_1 table.
    !> Expects w in primitive form. Uses FI bypass for hot cells.
    subroutine mhd_get_csound2_lte(w, x, ixI^L, ixO^L, cs2)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: cs2(ixi^s)
 
        double precision :: nh_val, log_nh, log_p_nh, g1, p_over_rho
        integer :: ix^d
 
        timeeos0 = mpi_wtime()
 
        if (eos%gamma1_method == 'constant') then
            cs2(ixo^s) = eos%gamma * w(ixo^s, p_) / w(ixo^s, rho_)
        else
            {do ix^db=ixomin^db,ixomax^db\}
                p_over_rho = w(ix^d, p_) / w(ix^d, rho_)
                if (p_over_rho > eos%p_rho_FI_threshold) then
                    cs2(ix^d) = eos%gamma * p_over_rho
                else
                    nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
                    if (eos%method == 'analytic') then
                        if (iw_te > 0 .and. w(ix^d,iw_te) > 0.0d0) then
                            g1 = saha_gamma1_from_nh_t(nh_val, w(ix^d,iw_te))
                        else
                            g1 = eos%gamma
                        end if
                    else
                        log_nh = dlog10(nh_val)
                        log_p_nh = dlog10(w(ix^d, p_) / nh_val)
                        g1 = gamma1_from_nh_p(log_nh, log_p_nh)
                    end if
                    cs2(ix^d) = g1 * p_over_rho
                end if
            {end do\}
        end if
 
        timeeos_csound = timeeos_csound + (mpi_wtime()-timeeos0)
 
    end subroutine mhd_get_csound2_lte
 
    subroutine mhd_get_gamma1_lte(w, x, ixI^L, ixO^L, gamma1)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(in)    :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
        double precision, intent(out)   :: gamma1(ixi^s)
 
        double precision :: nh_val, p_over_rho
        integer :: ix^d
 
        if (eos%gamma1_method == 'constant') then
            gamma1(ixo^s) = eos%gamma
            return
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            p_over_rho = w(ix^d, p_) / w(ix^d, rho_)
            if (p_over_rho > eos%p_rho_FI_threshold) then
                gamma1(ix^d) = eos%gamma
            else
                nh_val = w(ix^d, rho_) / eos%nH2rhoFactor
                if (eos%method == 'analytic') then
                    if (iw_te > 0 .and. w(ix^d,iw_te) > 0.0d0) then
                        gamma1(ix^d) = saha_gamma1_from_nh_t(nh_val, w(ix^d,iw_te))
                    else
                        gamma1(ix^d) = eos%gamma
                    end if
                else
                    gamma1(ix^d) = gamma1_from_nh_p(dlog10(nh_val), &
                        dlog10(w(ix^d, p_) / nh_val))
                end if
            end if
        {end do\}
 
    end subroutine mhd_get_gamma1_lte
 
    !> Rfactor routines (matching HD: EoS functions live in EoS module)
 
    !> Rfactor = p/(rho*T) for constant ionisation degree (FI/PI no-energy).
    !> Stays in the seam: RR is the physics module's gas-constant factor, not
    !> visible to mod_eos. Rfactor_from_LTE lives in mod_eos (no RR).
    subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
        use mod_global_parameters
        integer, intent(in) :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,1:nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: rfactor(ixi^s)
 
        rfactor(ixo^s)=rr
 
    end subroutine rfactor_from_constant_ionization
 
    !> Rfactor from partial ionisation temperature lookup (Leenaarts et al. 2012)
 
    !> Thermal pressure routines (moved from mod_mhd_phys.t)
 
    !> Calculate isothermal thermal pressure
    subroutine mhd_get_pthermal_noe(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        if(has_equi_rho_and_p) then
            pth(ixo^s)=mhd_adiab*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,0))**eos%gamma
        else
            pth(ixo^s)=mhd_adiab*w(ixo^s,rho_)**eos%gamma
        end if
 
    end subroutine mhd_get_pthermal_noe
 
    !> Calculate thermal pressure from internal energy
    subroutine mhd_get_pthermal_inte(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        use mod_small_values, only: trace_small_values
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        integer :: iw, ix^d
 
        {do ix^db= ixomin^db,ixomax^db\}
            if(has_equi_rho_and_p) then
                pth(ix^d)=eos%gamma_minus_1*w(ix^d,e_)+block%equi_vars(ix^d,equi_pe0_,0)
            else
                pth(ix^d)=eos%gamma_minus_1*w(ix^d,e_)
            end if
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
 
        if(check_small_values.and..not.fix_small_values) then
            {do ix^db= ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                         " encountered when call mhd_get_pthermal_inte"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {end do\}
        end if
 
    end subroutine mhd_get_pthermal_inte
 
    !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho-b**2/2) within ixO^L
    subroutine mhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        use mod_small_values, only: trace_small_values
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        integer :: iw, ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            if(has_equi_rho_and_p) then
                pth(ix^d)=eos%gamma_minus_1*(w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0))&
                     +(^c&w(ix^d,b^c_)**2+)))+block%equi_vars(ix^d,equi_pe0_,0)
            else
                pth(ix^d)=eos%gamma_minus_1*(w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_)&
                     +(^c&w(ix^d,b^c_)**2+))-mhd_uawsom_wave_energy_cell(w(ix^d,:)))
            end if
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
 
        if(check_small_values.and..not.fix_small_values) then
            {do ix^db=ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                         " encountered when call mhd_get_pthermal"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {end do\}
        end if
 
    end subroutine mhd_get_pthermal_origin
 
    !> Calculate thermal pressure for LTE EoS (delegates to eos%get_thermal_pressure)
    subroutine mhd_get_pthermal_lte(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        use mod_small_values, only: trace_small_values
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
        double precision :: nh(ixi^s)
 
        integer :: iw, ix^d
 
        ! LTE: p = nH * (1 + He + ne/nH) * T from stored state variables
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        pth(ixo^s) = nh(ixo^s) * (1.0d0 + eos%He_abundance &
            + (w(ixo^s,ne_) / nh(ixo^s))) * w(ixo^s,te_)
 
        if(fix_small_values) then
            {do ix^db=ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
            {end do\}
        else if(check_small_values) then
            {do ix^db=ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                         " encountered when call mhd_get_pthermal_LTE"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {end do\}
        end if
 
    end subroutine mhd_get_pthermal_lte
 
    !> Calculate thermal pressure for semirelativistic MHD
    subroutine mhd_get_pthermal_semirelati(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        use mod_small_values, only: trace_small_values
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        double precision :: b(ixo^s,1:ndir), v(ixo^s,1:ndir), tmp, b2, gamma2, inv_rho
        integer :: iw, ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            b2=(^c&w(ix^d,b^c_)**2+)
            if(b2>smalldouble) then
                tmp=1.d0/sqrt(b2)
            else
                tmp=0.d0
            end if
            ^c&b(ix^d,^c)=w(ix^d,b^c_)*tmp\
            tmp=(^c&b(ix^d,^c)*w(ix^d,m^c_)+)
 
            inv_rho=1.d0/w(ix^d,rho_)
            ! Va^2/c^2
            b2=b2*inv_rho*eos%inv_squared_c
            ! equation (15)
            gamma2=1.d0/(1.d0+b2)
            ! Convert momentum to velocity
            ^c&v(ix^d,^c)=gamma2*(w(ix^d,m^c_)+b2*b(ix^d,^c)*tmp)*inv_rho\
 
            ! E=Bxv
            {^ifthreec
                b(ix^d,1)=w(ix^d,b2_)*v(ix^d,3)-w(ix^d,b3_)*v(ix^d,2)
                b(ix^d,2)=w(ix^d,b3_)*v(ix^d,1)-w(ix^d,b1_)*v(ix^d,3)
                b(ix^d,3)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
            }
            {^iftwoc
                b(ix^d,1)=zero
                b(ix^d,2)=w(ix^d,b1_)*v(ix^d,2)-w(ix^d,b2_)*v(ix^d,1)
            }
            {^ifonec
                b(ix^d,1)=zero
            }
            ! Calculate pressure = (gamma-1) * (e-eK-eB-eE)
            pth(ix^d)=eos%gamma_minus_1*(w(ix^d,e_)&
                       -half*((^c&v(ix^d,^c)**2+)*w(ix^d,rho_)&
                       +(^c&w(ix^d,b^c_)**2+)&
                       +(^c&b(ix^d,^c)**2+)*eos%inv_squared_c))
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
 
        if(check_small_values.and..not.fix_small_values) then
            {do ix^db=ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                         " encountered when call mhd_get_pthermal_semirelati"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {end do\}
        end if
 
    end subroutine mhd_get_pthermal_semirelati
 
    !> Calculate thermal pressure=(gamma-1)*(e-0.5*m**2/rho) within ixO^L
    subroutine mhd_get_pthermal_hde(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        use mod_small_values, only: trace_small_values
 
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        integer :: iw, ix^d
 
        {do ix^db= ixomin^db,ixomax^db\}
            pth(ix^d)=eos%gamma_minus_1*(w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_)))
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
        if(check_small_values.and..not.fix_small_values) then
            {do ix^db= ixomin^db,ixomax^db\}
                if(pth(ix^d)<small_pressure) then
                    write(*,*) "Error: small value of gas pressure",pth(ix^d),&
                         " encountered when call mhd_get_pthermal_hde"
                    write(*,*) "Iteration: ", it, " Time: ", global_time
                    write(*,*) "Location: ", x(ix^d,:)
                    write(*,*) "Cell number: ", ix^d
                    do iw=1,nw
                        write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
                    end do
                    if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
                    write(*,*) "Saving status at the previous time step"
                    crash=.true.
                end if
            {end do\}
        end if
 
    end subroutine mhd_get_pthermal_hde
 
    !> Temperature routines (moved from mod_mhd_phys.t)
 
    !> Copy temperature from stored Te variable
    subroutine mhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
        res(ixo^s) = w(ixo^s, te_)
    end subroutine mhd_get_temperature_from_te
 
    !> Calculate temperature=p/rho when in e_ the internal energy is stored
    subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: r(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,r)
        res(ixo^s) = eos%gamma_minus_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
    end subroutine mhd_get_temperature_from_eint
 
    !> Calculate temperature=p/rho when in e_ the total energy is stored
    subroutine mhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: r(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,r)
        call mhd_get_pthermal(w,x,ixi^l,ixo^l,res)
        res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
 
    end subroutine mhd_get_temperature_from_etot
 
    subroutine mhd_get_temperature_from_etot_with_equi(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: r(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,r)
        call mhd_get_pthermal(w,x,ixi^l,ixo^l,res)
        res(ixo^s)=res(ixo^s)/(r(ixo^s)*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,b0i)))
 
    end subroutine mhd_get_temperature_from_etot_with_equi
 
    subroutine mhd_get_temperature_from_eint_with_equi(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: r(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,r)
        res(ixo^s) = (eos%gamma_minus_1 * w(ixo^s, e_) + block%equi_vars(ixo^s,equi_pe0_,b0i)) /&
                    ((w(ixo^s,rho_) +block%equi_vars(ixo^s,equi_rho0_,b0i))*r(ixo^s))
 
    end subroutine mhd_get_temperature_from_eint_with_equi
 
    !> LTE+MHD: T from total energy. Subtract KE+ME (via the MHD-aware
    !> phys_e_to_ei dispatcher) to get eint, then look up T in the LTE table
    !> through eos%get_temperature_from_eint (bound to get_temperature_from_eint_LTE
    !> by eos_finalise). Always fresh; does not use the iw_te cache.
    !> Should NOT be bound when mhd_internal_e=.true. (w(:,e_) is already eint
    !> and phys_e_to_ei is unbound) - use eos%get_temperature_from_eint directly.
    subroutine mhd_get_temperature_from_etot_lte(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: wlocal(ixi^s, 1:nw)
 
        wlocal(ixi^s, 1:nw) = w(ixi^s, 1:nw)
        ! Subtract KE+ME from w(:,e_) via the MHD-aware dispatcher
        call phys_e_to_ei(ixi^l, ixo^l, wlocal, x)
        ! Now wlocal(:,e_) is the gas internal energy; LTE table lookup
        call eos%get_temperature_from_eint(wlocal, x, ixi^l, ixo^l, res)
 
    end subroutine mhd_get_temperature_from_etot_lte
 
    subroutine mhd_get_temperature_equi(w,x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
 
        double precision :: r(ixi^s)
 
        call eos%get_Rfactor(w,x,ixi^l,ixo^l,r)
        res(ixo^s)= block%equi_vars(ixo^s,equi_pe0_,b0i)/(block%equi_vars(ixo^s,equi_rho0_,b0i)*r(ixo^s))
 
    end subroutine mhd_get_temperature_equi
 
    !> Equilibrium extraction routines (moved from mod_mhd_phys.t)
 
    subroutine mhd_get_rho_equi(w, x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
        res(ixo^s) = block%equi_vars(ixo^s,equi_rho0_,b0i)
    end subroutine mhd_get_rho_equi
 
    subroutine mhd_get_pe_equi(w,x, ixI^L, ixO^L, res)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, 1:nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: res(ixi^s)
        res(ixo^s) = block%equi_vars(ixo^s,equi_pe0_,b0i)
    end subroutine mhd_get_pe_equi
 
    !> Electron and hydrogen number densities of the background alone, for the
    !> equilibrium rate that subtract_equi removes. Splitting is restricted to
    !> the FI EoS, so the electron count is the fully ionised one.
    subroutine mhd_get_ne_nh_equi(ixI^L, ixO^L, w, x, ne, nH)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, nw)
        double precision, intent(in) :: x(ixi^s, 1:ndim)
        double precision, intent(out):: ne(ixi^s), nh(ixi^s)
        nh(ixo^s) = block%equi_vars(ixo^s,equi_rho0_,b0i) / eos%nH2rhoFactor
        ne(ixo^s) = nh(ixo^s) * eos%neOnH_FI
    end subroutine mhd_get_ne_nh_equi
 
    !> Internal energy extraction + small value handling (moved from mod_mhd_phys.t)
 
    !> Get internal energy from conserved state (origin formulation)
    function mhd_get_ei_origin(w, ixI^L, ixO^L) result(ei)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, nw)
        double precision             :: ei(ixo^s)
 
        ! ei = E_total - e_kinetic - e_magnetic
        ei(ixo^s) = w(ixo^s,e_) - half*((^c&w(ixo^s,m^c_)**2+)/w(ixo^s,rho_) &
             + (^c&w(ixo^s,b^c_)**2+))
    end function mhd_get_ei_origin
 
    !> Get internal energy from conserved state (internal energy formulation)
    function mhd_get_ei_inte(w, ixI^L, ixO^L) result(ei)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s, nw)
        double precision             :: ei(ixo^s)
 
        ! e_ already stores internal energy
        ei(ixo^s) = w(ixo^s,e_)
    end function mhd_get_ei_inte
 
    ! mhd_handle_small_ei stays in mod_mhd_phys.t (callers are there)
 
    !> Update temperature Te_ from pressure using partial ionization
 
    !> PI energy-mode routines (eos_type='PI', ionE=.true.)
    !> Structural mirror of the FI conversion/pthermal/csound family, sharing
    !> FI's eq_state_units / RR=1 normalisation; the only difference is that the
    !> eint<->p relation is delegated to the portable scalar backend
    !> (mod_eos_PI), so eint carries the ionisation-energy term.
    !> Origin: w(e_) is total energy; inte: w(e_) is gas internal energy.
 
    !> Primitive pressure -> total energy (origin). m^C_ is still velocity here
    !> (momentum conversion follows in to_conserved), so KE = half*rho*v^2.
    subroutine mhd_p_to_e_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: eint
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            call eint_from_rho_p_pi(w(ix^d,rho_), w(ix^d,p_), eint)
            w(ix^d,e_)=eint &
                      +half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                      +(^c&w(ix^d,b^c_)**2+))
        {end do\}
 
    end subroutine mhd_p_to_e_pi
 
    !> Primitive pressure -> gas internal energy (inte formulation)
    subroutine mhd_p_to_eint_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: eint
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            call eint_from_rho_p_pi(w(ix^d,rho_), w(ix^d,p_), eint)
            w(ix^d,e_)=eint
        {end do\}
 
    end subroutine mhd_p_to_eint_pi
 
    !> Primitive -> conserved (origin, total energy)
    subroutine mhd_to_conserved_origin_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        call mhd_p_to_e_pi(ixi^l, ixo^l, w, x)
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
 
    end subroutine mhd_to_conserved_origin_pi
 
    !> Primitive -> conserved (internal energy)
    subroutine mhd_to_conserved_inte_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        call mhd_p_to_eint_pi(ixi^l, ixo^l, w, x)
        ^c&w(ixo^s,m^c_)=w(ixo^s,rho_)*w(ixo^s,m^c_)\
 
    end subroutine mhd_to_conserved_inte_pi
 
    !> Conserved -> primitive (origin). Subtract KE+ME -> eint, invert -> p.
    subroutine mhd_to_primitive_origin_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, t, rfac
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'mhd_to_primitive_origin_PI')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ! Convert momentum to velocity
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            ! eint = E - 0.5*rho*v^2 - 0.5*B^2
            eint_val = w(ix^d,e_) &
                      - half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+) &
                      + (^c&w(ix^d,b^c_)**2+))
            eint_val = max(eint_val, smalldouble)
            call state_from_eint_pi(w(ix^d,rho_), eint_val, t, w(ix^d,p_), rfac)
        {end do\}
 
    end subroutine mhd_to_primitive_origin_pi
 
    !> Conserved -> primitive (internal energy). e_ already holds eint.
    subroutine mhd_to_primitive_inte_pi(ixI^L,ixO^L,w,x)
        use mod_global_parameters
        integer, intent(in)             :: ixi^l, ixo^l
        double precision, intent(inout) :: w(ixi^s, nw)
        double precision, intent(in)    :: x(ixi^s, 1:ndim)
 
        double precision :: inv_rho, eint_val, t, rfac
        integer :: ix^d
 
        if (fix_small_values) then
            call phys_handle_small_values(.false., w, x, ixi^l, ixo^l, &
                'mhd_to_primitive_inte_PI')
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            eint_val = max(w(ix^d,e_), smalldouble)
            call state_from_eint_pi(w(ix^d,rho_), eint_val, t, w(ix^d,p_), rfac)
            ! Convert momentum to velocity
            inv_rho = 1.d0/w(ix^d,rho_)
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
        {end do\}
 
    end subroutine mhd_to_primitive_inte_pi
 
    !> Thermal pressure (origin): subtract KE+ME, invert eint -> p.
    subroutine mhd_get_pthermal_origin_pi(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        double precision :: eint_val, t, rfac
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            eint_val = w(ix^d,e_) &
                      - half*((^c&w(ix^d,m^c_)**2+)/w(ix^d,rho_) &
                      + (^c&w(ix^d,b^c_)**2+))
            eint_val = max(eint_val, smalldouble)
            call state_from_eint_pi(w(ix^d,rho_), eint_val, t, pth(ix^d), rfac)
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
 
    end subroutine mhd_get_pthermal_origin_pi
 
    !> Thermal pressure (inte): e_ holds eint, invert eint -> p.
    subroutine mhd_get_pthermal_inte_pi(w,x,ixI^L,ixO^L,pth)
        use mod_global_parameters
        integer, intent(in)          :: ixi^l, ixo^l
        double precision, intent(in) :: w(ixi^s,nw)
        double precision, intent(in) :: x(ixi^s,1:ndim)
        double precision, intent(out):: pth(ixi^s)
 
        double precision :: eint_val, t, rfac
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            eint_val = max(w(ix^d,e_), smalldouble)
            call state_from_eint_pi(w(ix^d,rho_), eint_val, t, pth(ix^d), rfac)
            if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
        {end do\}
 
    end subroutine mhd_get_pthermal_inte_pi
 
    !> Adiabatic sound speed squared from primitive (rho, p).
 
    !> Effective Gamma1 = cs2 * rho / p (for the same primitive state).
 
    !> PI energy mode: refresh Te_ from gas internal energy via the backend
    !> eint->T inversion (direct, so no lagged wCT -- unlike no-energy
    !> mhd_update_temperature which lags iz_H via wCT(Te_)). Uses MHD's Te_ index
    !> (PI registers Te via var_set_auxvar, so the generic iw_te is unset).
 
    !> PI prominence (T,p) routines (eos_type='PI', pi_table='prominence').
    !> Prominence ionisation depends on (T,p), not T alone, so the chromosphere
    !> T-only seam does not apply. Following the collaborator's design, R and T
    !> are recomputed from (rho, pth) via ionization_get_state, which dispatches
    !> the prominence (pressure-bisection) inversion internally. No-energy only
    !> (guarded in eos_validate_params), so pth = (gamma-1)*eint is R-independent
    !> and obtained from phys_get_ei (conserved w; KE/ME removed).
 
    !> R-factor for the prominence (T,p) table: recompute from (rho, pth).
 
    !> Refresh stored Te_ for the prominence table from (rho, pth).
 
end module mod_mhd_eos
!> Needs a line after to pass the preprocesor