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_radiative_cooling, only: build_y_mod_table, rc_t_ymod_max_k
 
    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
                if(eos%ionE) call mpistop('PI energy EoS not supported with mhd_hydrodynamic_e')
                eos%to_primitive        => mhd_to_primitive_hde
                eos%to_conserved        => mhd_to_conserved_hde
                phys_get_pthermal        => mhd_get_pthermal_hde
                mhd_get_pthermal         => mhd_get_pthermal_hde
            else if(mhd_semirelativistic) then
                if(mhd_energy) then
                    ! the eint <-> p step inside dispatches on the EoS
                    eos%to_primitive        => mhd_to_primitive_semirelati
                    eos%to_conserved        => mhd_to_conserved_semirelati
                    if(mhd_internal_e .and. eos%ionE) then
                      phys_get_pthermal => mhd_get_pthermal_inte_pi
                    else if(mhd_internal_e) then
                      phys_get_pthermal => mhd_get_pthermal_inte
                    else
                      phys_get_pthermal => mhd_get_pthermal_semirelati
                    end if
                    mhd_get_pthermal => phys_get_pthermal
                else
                    eos%to_primitive        => mhd_to_primitive_semirelati_noe
                    eos%to_conserved        => mhd_to_conserved_semirelati_noe
                    phys_get_pthermal       => mhd_get_pthermal_noe
                    mhd_get_pthermal        => mhd_get_pthermal_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
                    if(eos%ionE) 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_primitive  => mhd_to_primitive_inte
                      eos%to_conserved  => mhd_to_conserved_inte
                      phys_get_pthermal => mhd_get_pthermal_inte
                      mhd_get_pthermal  => mhd_get_pthermal_inte
                    end if
                else if(mhd_energy) then
                    if(eos%ionE) then
                      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
                    else
                      eos%to_primitive  => mhd_to_primitive_origin
                      eos%to_conserved  => mhd_to_conserved_origin
                      phys_get_pthermal => mhd_get_pthermal_origin
                      mhd_get_pthermal  => mhd_get_pthermal_origin
                    end if
                else
                    eos%to_primitive    => mhd_to_primitive_origin_noe
                    eos%to_conserved    => mhd_to_conserved_origin_noe
                    phys_get_pthermal   => mhd_get_pthermal_noe
                    mhd_get_pthermal    => mhd_get_pthermal_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 (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_semirelativistic) then
                eos%to_primitive  => mhd_to_primitive_semirelati
                eos%to_conserved  => mhd_to_conserved_semirelati
            else 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 if (eos%eos_type == 'PI' .and. eos%ionE) then
            eos%get_csound2 => mhd_get_csound2_fi
            phys_get_gamma1 => get_gamma1_pi
        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
 
        ! 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
            if(mhd_semirelativistic) then
              phys_get_ei => mhd_get_ei_semirelati
            else
              phys_get_ei => mhd_get_ei_origin
            end if
        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
        phys_to_primitive        => eos%to_primitive
        phys_to_conserved        => eos%to_conserved
        eos%get_thermal_pressure => mhd_get_pthermal
 
    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)
        if (eos%eos_type == 'LTE' .or. (eos%eos_type == 'PI' .and. eos%ionE)) then
            ! eint carries the ionisation energy, so T comes from the EoS's own
            ! eint->T inversion (LTE tables, or the PI energy-mode backend)
            if (mhd_internal_e) then
                ! w(:,e_) is eint; both hooks do the EoS lookup directly
                eos%get_temperature_from_etot => eos%get_temperature_from_eint
            else
                ! Total-energy: subtract KE+ME (via phys_e_to_ei), then the EoS eint->T lookup
                eos%get_temperature_from_etot => mhd_get_temperature_from_etot_lte
            end if
        else  ! FI, or PI without ionisation energy
            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
            call eos_bind_tc(tc_fl, .not. associated(eos%get_Rfactor, rfactor_from_constant_ionization))
            ! 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_analytic   =  (eos%method_id == eos_analytic)
            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%log_T_floor       =  eos_get_log_t_floor()
            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
            rc_fl%T_Ymod_max        =  rc_t_ymod_max_k/unit_temperature
            ! 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, any EoS)
    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 :: efd(1:ndir), s(1:ndir), eint(ixo^s)
        double precision :: btot(1:ndir)
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! this loop can not be vectorized
            eint(ix^d)=mhd_semirelati_eint(w(ix^d,p_),w(ix^d,rho_))
        {end do\}
        {do ix^db=ixomin^db,ixomax^db\}
            ! fully vectorized loop, E, S, e_, m^C_ all in one loop for speed
            ^c&btot(^c)=w(ix^d,b^c_)\
            if(b0field) then
                ^c&btot(^c)=btot(^c)+block%B0(ix^d,^c,b0i)\
            end if
            {^ifthreec
                efd(1)=btot(2)*w(ix^d,m3_)-btot(3)*w(ix^d,m2_)
                efd(2)=btot(3)*w(ix^d,m1_)-btot(1)*w(ix^d,m3_)
                efd(3)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                s(1)=efd(2)*btot(3)-efd(3)*btot(2)
                s(2)=efd(3)*btot(1)-efd(1)*btot(3)
                s(3)=efd(1)*btot(2)-efd(2)*btot(1)
            }
            {^iftwoc
                efd(1)=zero
                ! switch 3 with 2 to add 3 when ^C from 1 to 2
                efd(2)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                s(1)=-efd(2)*btot(2)
                s(2)=efd(2)*btot(1)
            }
            {^ifonec
                efd(1)=zero
                s(1)=zero
            }
            if(mhd_internal_e) then
                w(ix^d,e_)=eint(ix^d)
            else
                ! equation (9), the magnetic energy from the stored slot (eq 92 under splitting)
                w(ix^d,e_)=eint(ix^d)+half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                          +(^c&w(ix^d,b^c_)**2+)&
                          +(^c&efd(^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(^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 :: efd(1:ndir), s(1:ndir)
        double precision :: btot(1:ndir)
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! fully vectorized loop, E, S, e_, m^C_ all in one loop for speed
            ^c&btot(^c)=w(ix^d,b^c_)\
            if(b0field) then
                ^c&btot(^c)=btot(^c)+block%B0(ix^d,^c,b0i)\
            end if
            {^ifthreec
                efd(1)=btot(2)*w(ix^d,m3_)-btot(3)*w(ix^d,m2_)
                efd(2)=btot(3)*w(ix^d,m1_)-btot(1)*w(ix^d,m3_)
                efd(3)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                s(1)=efd(2)*btot(3)-efd(3)*btot(2)
                s(2)=efd(3)*btot(1)-efd(1)*btot(3)
                s(3)=efd(1)*btot(2)-efd(2)*btot(1)
            }
            {^iftwoc
                efd(1)=zero
                ! switch 3 with 2 to add 3 when ^C from 1 to 2
                efd(2)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                s(1)=-efd(2)*btot(2)
                s(2)=efd(2)*btot(1)
            }
            {^ifonec
                efd(1)=zero
                s(1)=zero
            }
            ! Convert velocity to momentum, equation (9)
            ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)+s(^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, any EoS)
    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 :: btot(1:ndir), efd(1:ndir)
        double precision :: b2, tmp, inv_rho, rhoc2
        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\}
            ^c&btot(^c)=w(ix^d,b^c_)\
            if(b0field) then
              ^c&btot(^c)=btot(^c)+block%B0(ix^d,^c,b0i)\
            end if
            ! Keppens's formula to convert momentum to velocity
            b2=(^c&btot(^c)**2+)
            tmp=(^c&btot(^c)*w(ix^d,m^c_)+)
            rhoc2=w(ix^d,rho_)*mhd_squared_c
            inv_rho=1.d0/(w(ix^d,rho_)*(rhoc2+b2))
            ^c&w(ix^d,m^c_)=(rhoc2*w(ix^d,m^c_)+btot(^c)*tmp)*inv_rho\
            if(.not.mhd_internal_e) then
              ! E=Bxv
              {^ifthreec
              efd(1)=btot(2)*w(ix^d,m3_)-btot(3)*w(ix^d,m2_)
              efd(2)=btot(3)*w(ix^d,m1_)-btot(1)*w(ix^d,m3_)
              efd(3)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
              }
              {^iftwoc
              efd(1)=zero
              efd(2)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
              }
              {^ifonec
              efd(1)=zero
              }
              w(ix^d,e_)=w(ix^d,e_)-half*((^c&w(ix^d,m^c_)**2+)*w(ix^d,rho_)&
                         +(^c&w(ix^d,b^c_)**2+)&
                         +(^c&efd(^c)**2+)*eos%inv_squared_c)
            end if
        {end do\}
        {do ix^db=ixomin^db,ixomax^db\}
            ! loop can not be vectorized
            w(ix^d,p_)=mhd_semirelati_p(w(ix^d,e_),w(ix^d,rho_))
        {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 :: btot(1:ndir)
        double precision :: b2, rhoc2, tmp, 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\}
            ^c&btot(^c)=w(ix^d,b^c_)\
            if(b0field) then
              ^c&btot(^c)=btot(^c)+block%B0(ix^d,^c,b0i)\
            end if
            ! Keppens's formula to convert momentum to velocity
            b2=(^c&btot(^c)**2+)
            tmp=(^c&btot(^c)*w(ix^d,m^c_)+)
            rhoc2=w(ix^d,rho_)*mhd_squared_c
            inv_rho=1.d0/(w(ix^d,rho_)*(rhoc2+b2))
            ^c&w(ix^d,m^c_)=(rhoc2*w(ix^d,m^c_)+btot(^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)
 
        ! 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_)\
 
    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)
 
        ! 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_)\
 
    end subroutine mhd_to_conserved_inte_lte
 
    !> Convert pressure to total energy (eint + eK + eB) for origin formulation
    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
 
        call mhd_p_to_eint(ixi^l,ixo^l,w,x)
        {do ix^db=ixomin^db,ixomax^db\}
            w(ix^d,e_) = w(ix^d,e_) + 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_)
        {end do\}
 
    end subroutine mhd_p_to_e
 
    !> Convert pressure to internal energy only (for internal_e formulation)
    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 :: nh(ixi^s), log_nh(ixi^s)
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        log_nh(ixo^s) = 0.0d0
        if (eos%ionE) log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            w(ix^d,e_) = lte_eint_from_p(w(ix^d,p_), nh(ix^d), log_nh(ix^d), w(ix^d,rho_))
        {end do\}
 
    end subroutine mhd_p_to_eint
 
    !> LTE: Conserved (rho, rho*v, E_total, B) -> Primitive (rho, v, p, B)
    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
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        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)
        log_nh(ixo^s) = 0.0d0
        if (eos%ionE) log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
            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+))
            w(ix^d,p_) = lte_p_from_eint(eint_val, nh(ix^d), log_nh(ix^d), inv_rho)
        {end do\}
 
    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
        double precision :: nh(ixi^s), log_nh(ixi^s)
        integer :: ix^d
 
        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)
        log_nh(ixo^s) = 0.0d0
        if (eos%ionE) log_nh(ixo^s) = dlog10(nh(ixo^s))
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            w(ix^d,p_) = lte_p_from_eint(w(ix^d,e_), nh(ix^d), log_nh(ix^d), inv_rho)
            ^c&w(ix^d,m^c_)=w(ix^d,m^c_)*inv_rho\
        {end do\}
 
    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)
        double precision :: vel(ixo^s,1:ndir), emech(ixo^s)
        integer :: ix^d
 
        call eos%get_nH(w, x, ixi^l, ixo^l, nh)
        log_nh(ixo^s) = dlog10(nh(ixo^s))
        if (mhd_semirelativistic) call mhd_semirelati_vel_emech(w,ixi^l,ixo^l,vel,emech)
 
        {do ix^db=ixomin^db,ixomax^db\}
            inv_rho = 1.d0/w(ix^d,rho_)
            if (mhd_semirelativistic) then
                ! momentum and energy carry the electric-field terms
                ^c&w(ix^d,m^c_)=vel(ix^d,^c)\
                eint_val = w(ix^d,e_) - emech(ix^d)
            else
                ! 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
            end if
            eint_val = max(eint_val, nh(ix^d)*eos%eint_nH_min, 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_id == eos_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(ixo^s), t_fi, log_t_min
        double precision :: nh(ixi^s), btot(1:ndir)
        double precision :: efd(1:ndir), s(1:ndir)
        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_id == eos_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)
 
        {do ix^db=ixomin^db,ixomax^db\}
            ! loop can not be vectorized
            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(ix^d) = nh(ix^d) &
                    * (eos%n_per_nH_FI * t_val * eos%inv_gamma_minus_1 &
                    + eos%eion_per_nH)
            else if (eos%method_id == eos_analytic) then
                ! Analytical Saha: eint from T directly
                eint_val(ix^d) = saha_eint_from_nh_t(nh(ix^d), t_val) * nh(ix^d)
            else
                ! Ionisation zone: eint/nH from T table
                eint_val(ix^d) = eint_nh_from_t( &
                    dlog10(nh(ix^d)), &
                    dlog10(max(t_val, 10.0d0**log_t_min))) &
                    * nh(ix^d)
            end if
        {end do\}
 
        if (mhd_semirelativistic) then
          {do ix^db=ixomin^db,ixomax^db\}
            ! fully vectorized loop, E, S, e_, m^C_ all in one loop for speed
              ^c&btot(^c)=w(ix^d,b^c_)\
              if(b0field) then
                  ^c&btot(^c)=btot(^c)+block%B0(ix^d,^c,b0i)\
              end if
              {^ifthreec
                  efd(1)=btot(2)*w(ix^d,m3_)-btot(3)*w(ix^d,m2_)
                  efd(2)=btot(3)*w(ix^d,m1_)-btot(1)*w(ix^d,m3_)
                  efd(3)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                  s(1)=efd(2)*btot(3)-efd(3)*btot(2)
                  s(2)=efd(3)*btot(1)-efd(1)*btot(3)
                  s(3)=efd(1)*btot(2)-efd(2)*btot(1)
              }
              {^iftwoc
                  efd(1)=zero
                  ! switch 3 with 2 to add 3 when ^C from 1 to 2
                  efd(2)=btot(1)*w(ix^d,m2_)-btot(2)*w(ix^d,m1_)
                  s(1)=-efd(2)*btot(2)
                  s(2)=efd(2)*btot(1)
              }
              {^ifonec
                  efd(1)=zero
                  s(1)=zero
              }
              if (mhd_internal_e) then
                  w(ix^d,e_) = eint_val(ix^d)
              else
                  ! equation (9), the magnetic energy from the stored slot
                  w(ix^d,e_) = eint_val(ix^d) &
                      + half*(w(ix^d,rho_)*(^c&w(ix^d,m^c_)**2+) &
                      + (^c&w(ix^d,b^c_)**2+) &
                      + (^c&efd(^c)**2+)*eos%inv_squared_c)
              end if
              ! Convert velocity to momentum
              ^c&w(ix^d,m^c_)=w(ix^d,rho_)*w(ix^d,m^c_)+s(^c)*eos%inv_squared_c\
          {end do\}
        else
          {do ix^db=ixomin^db,ixomax^db\}
              if (mhd_internal_e) then
                  w(ix^d,e_) = eint_val(ix^d)
              else
                  ! E = eint + 0.5*rho*v^2 + 0.5*B^2
                  w(ix^d,e_) = eint_val(ix^d) &
                      + 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 if
 
    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, ploc
        integer :: ix^d
 
        {do ix^db=ixomin^db,ixomax^db \}
            if(has_equi_rho_and_p) then
              rho=(w(ix^d,rho_)+block%equi_vars(ix^d,equi_rho0_,0))
              ploc=(w(ix^d,p_)+block%equi_vars(ix^d,equi_pe0_,0))
              ! sound speed**2 
              cs2(ix^d)=eos%gamma*ploc/rho
            else
              ! sound speed**2 
              cs2(ix^d)=eos%gamma*w(ix^d,p_)/w(ix^d,rho_)
            end if
        {end do\}
 
    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)
 
        integer :: ix^d
 
        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\}
                cs2(ix^d) = lte_gamma1_from_p(w(ix^d,p_), w(ix^d,rho_)) * (w(ix^d,p_) / w(ix^d,rho_))
            {end do\}
        end if
 
    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)
 
        integer :: ix^d
 
        if (eos%gamma1_method == 'constant') then
            gamma1(ixo^s) = eos%gamma
            return
        end if
 
        {do ix^db=ixomin^db,ixomax^db\}
            gamma1(ix^d) = lte_gamma1_from_p(w(ix^d,p_), w(ix^d,rho_))
        {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
        use mod_usr_methods, only: usr_set_adiab, usr_set_gamma
 
        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 :: adiabs(ixi^s), gammas(ixi^s)
 
        if(associated(usr_set_adiab)) then
           call usr_set_adiab(w,x,ixi^l,ixo^l,adiabs)
        else
           adiabs=mhd_adiab
        end if
        if(associated(usr_set_gamma)) then
           call usr_set_gamma(w,x,ixi^l,ixo^l,gammas)
        else
           gammas=eos%gamma
        end if
 
        if(has_equi_rho_and_p) then
            pth(ixo^s)=adiabs(ixo^s)*(w(ixo^s,rho_)+block%equi_vars(ixo^s,equi_rho0_,0))**gammas(ixo^s)
        else
            pth(ixo^s)=adiabs(ixo^s)*w(ixo^s,rho_)**gammas(ixo^s)
        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 (any EoS)
    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 :: vel(ixo^s,1:ndir), emech(ixo^s)
        integer :: iw, ix^d
 
        if(mhd_internal_e) then
          {do ix^db=ixomin^db,ixomax^db\}
              pth(ix^d)=mhd_semirelati_p(w(ix^d,e_),w(ix^d,rho_))
              if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
          {end do\}
        else
          call mhd_semirelati_vel_emech(w,ixi^l,ixo^l,vel,emech)
          {do ix^db=ixomin^db,ixomax^db\}
              pth(ix^d)=mhd_semirelati_p(w(ix^d,e_)-emech(ix^d),w(ix^d,rho_))
              if(fix_small_values.and.pth(ix^d)<small_pressure) pth(ix^d)=small_pressure
          {end do\}
        end if
 
        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
 
    !> T from total energy for LTE, and for PI with ionisation energy. KE+ME are subtracted through
    !> the MHD-aware phys_e_to_ei, and T comes from the EoS eint->T lookup through
    !> eos%get_temperature_from_eint. Always fresh; does not use the iw_te cache. With
    !> mhd_internal_e w(:,e_) is already eint, and eos%get_temperature_from_eint is bound instead
    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 (semirelativistic formulation)
    function mhd_get_ei_semirelati(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)
 
        double precision :: vel(ixo^s,1:ndir), emech(ixo^s)
 
        call mhd_semirelati_vel_emech(w,ixi^l,ixo^l,vel,emech)
        ei(ixo^s)=w(ixo^s,e_)-emech(ixo^s)
    end function mhd_get_ei_semirelati
 
    !> 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