diff --git a/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl
new file mode 100644
index 00000000000..40f45c0e49b
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_mhd_shockcapturing_subcell.jl
@@ -0,0 +1,124 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations3D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+
+Weak magnetic blast wave setup taken from Section 6.1 of the paper:
+- A. M. Rueda-Ramírez, S. Hennemann, F. J. Hindenlang, A. R. Winters, G. J. Gassner (2021)
+  An entropy stable nodal discontinuous Galerkin method for the resistive MHD
+  equations. Part II: Subcell finite volume shock capturing
+  [doi: 10.1016/j.jcp.2021.110580](https://doi.org/10.1016/j.jcp.2021.110580)
+"""
+function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations3D)
+    # Center of the blast wave is selected for the domain [0, 3]^3
+    inicenter = (1.5, 1.5, 1.5)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    delta_0 = 0.1
+    r_0 = 0.3
+    lambda = exp(5.0 / delta_0 * (r - r_0))
+
+    p_inner = 0.9
+    p_outer = 5e-2
+    prim_inner = SVector(1.2, 0.1, 0.0, 0.1, p_inner, 1.0, 1.0, 1.0, 0.0)
+    prim_outer = SVector(1.2, 0.2, -0.4, 0.2, p_outer, 1.0, 1.0, 1.0, 0.0)
+    prim_vars = (prim_inner + lambda * prim_outer) / (1.0 + lambda)
+
+    return prim2cons(prim_vars, equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell_local_symmetric)
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell_local_symmetric)
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but with slightly less warping.
+# The mapping will be interpolated at tree level, and then refined without changing
+# the geometry interpolant.
+function mapping(xi_, eta_, zeta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+    zeta = 1.5 * zeta_ + 1.5
+
+    y = eta +
+        3 / 11 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+         cos(0.5 * pi * (2 * eta - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    x = xi +
+        3 / 11 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+         cos(2 * pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    z = zeta +
+        3 / 11 * (cos(0.5 * pi * (2 * x - 3) / 3) *
+         cos(pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    return SVector(x, y, z)
+end
+
+trees_per_dimension = (2, 2, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3,
+                 mapping = mapping,
+                 initial_refinement_level = 2,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
+
+cfl = 0.9
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+stage_callbacks = (SubcellLimiterIDPCorrection(),)
+
+sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
+                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  callback = callbacks);
diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl
index ee3dedd0545..3e53eb91634 100644
--- a/src/equations/ideal_glm_mhd_2d.jl
+++ b/src/equations/ideal_glm_mhd_2d.jl
@@ -1982,8 +1982,7 @@ end
 
 # State validation for Newton-bisection method of subcell IDP limiting
 @inline function Base.isvalid(u, equations::IdealGlmMhdEquations2D)
-    p = pressure(u, equations)
-    if u[1] <= 0 || p <= 0
+    if u[1] <= 0 || pressure(u, equations) <= 0
         return false
     end
     return true
diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl
index a433756a14b..620c3c1f2a3 100644
--- a/src/equations/ideal_glm_mhd_3d.jl
+++ b/src/equations/ideal_glm_mhd_3d.jl
@@ -272,7 +272,7 @@ Non-symmetric two-point flux discretizing the nonconservative (source) term of
 Powell and the Galilean nonconservative term associated with the GLM multiplier
 of the [`IdealGlmMhdEquations3D`](@ref).
 
-On curvilinear meshes, the implementation differs from the reference since we use the averaged 
+On curvilinear meshes, the implementation differs from the reference since we use the averaged
 normal direction for the metrics dealiasing.
 
 ## References
@@ -361,6 +361,189 @@ end
     return f
 end
 
+# For `VolumeIntegralSubcellLimiting` the nonconservative flux is created as a callable struct to
+# enable dispatch on the type of the nonconservative term (symmetric / jump).
+"""
+    flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
+                                                normal_direction::AbstractVector,
+                                                equations::IdealGlmMhdEquations3D)
+
+Non-symmetric two-point flux discretizing the nonconservative (source) term of
+Powell and the Galilean nonconservative term associated with the GLM multiplier
+of the [`IdealGlmMhdEquations3D`](@ref).
+
+This implementation uses a non-conservative term that can be written as the product
+of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm
+et al. [`flux_nonconservative_powell`](@ref) for conforming meshes but it yields different
+results on non-conforming meshes(!). On curvilinear meshes this formulation applies the
+local normal direction compared to the averaged one used in [`flux_nonconservative_powell`](@ref).
+
+The two other flux functions with the same name return either the local
+or symmetric portion of the non-conservative flux based on the type of the
+nonconservative_type argument, employing multiple dispatch. They are used to
+compute the subcell fluxes in dg_3d_subcell_limiters.jl.
+
+## References
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems.
+  [DOI: 10.1016/j.jcp.2023.112607](https://doi.org/10.1016/j.jcp.2023.112607).
+  [arXiv: 2211.14009](https://doi.org/10.48550/arXiv.2211.14009).
+"""
+@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr,
+                                                                         normal_direction::AbstractVector,
+                                                                         equations::IdealGlmMhdEquations3D)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
+
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+    # The factor 0.5 of the averages can be omitted since it is already applied when this
+    # function is called.
+    psi_avg = (psi_ll + psi_rr)
+    B1_avg = (B1_ll + B1_rr)
+    B2_avg = (B2_ll + B2_rr)
+    B3_avg = (B3_ll + B3_rr)
+
+    B_dot_n_avg = B1_avg * normal_direction[1] + B2_avg * normal_direction[2] +
+                  B3_avg * normal_direction[3]
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
+                 v3_ll * normal_direction[3]
+
+    # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+    # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2,3}, 0, 0, 0, v_{1,2,3})
+    f = SVector(0,
+                B1_ll * B_dot_n_avg,
+                B2_ll * B_dot_n_avg,
+                B3_ll * B_dot_n_avg,
+                v_dot_B_ll * B_dot_n_avg + v_dot_n_ll * psi_ll * psi_avg,
+                v1_ll * B_dot_n_avg,
+                v2_ll * B_dot_n_avg,
+                v3_ll * B_dot_n_avg,
+                v_dot_n_ll * psi_avg)
+
+    return f
+end
+
+"""
+    flux_nonconservative_powell_local_symmetric(u_ll, normal_direction_ll::AbstractVector,
+                                                equations::IdealGlmMhdEquations3D,
+                                                nonconservative_type::NonConservativeLocal,
+                                                nonconservative_term::Integer)
+
+Local part of the Powell and GLM non-conservative terms. Needed for the calculation of
+the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems. 
+  [DOI: 10.1016/j.jcp.2023.112607](https://doi.org/10.1016/j.jcp.2023.112607).
+  [arXiv: 2211.14009](https://doi.org/10.48550/arXiv.2211.14009).
+  
+This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl.
+"""
+@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll,
+                                                                         normal_direction_ll::AbstractVector,
+                                                                         equations::IdealGlmMhdEquations3D,
+                                                                         nonconservative_type::NonConservativeLocal,
+                                                                         nonconservative_term::Integer)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+
+    if nonconservative_term == 1
+        # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+        v1_ll = rho_v1_ll / rho_ll
+        v2_ll = rho_v2_ll / rho_ll
+        v3_ll = rho_v3_ll / rho_ll
+        v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+        f = SVector(0,
+                    B1_ll,
+                    B2_ll,
+                    B3_ll,
+                    v_dot_B_ll,
+                    v1_ll,
+                    v2_ll,
+                    v3_ll,
+                    0)
+    else # nonconservative_term == 2
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        v1_ll = rho_v1_ll / rho_ll
+        v2_ll = rho_v2_ll / rho_ll
+        v3_ll = rho_v3_ll / rho_ll
+        v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] +
+                     v3_ll * normal_direction_ll[3]
+
+        f = SVector(0,
+                    0,
+                    0,
+                    0,
+                    v_dot_n_ll * psi_ll,
+                    0,
+                    0,
+                    0,
+                    v_dot_n_ll)
+    end
+    return f
+end
+
+"""
+    flux_nonconservative_powell_local_symmetric(u_ll, normal_direction_avg::AbstractVector,
+                                                equations::IdealGlmMhdEquations3D,
+                                                nonconservative_type::NonConservativeSymmetric,
+                                                nonconservative_term::Integer)
+
+Symmetric part of the Powell and GLM non-conservative terms. Needed for the calculation of
+the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems.
+  [DOI: 10.1016/j.jcp.2023.112607](https://doi.org/10.1016/j.jcp.2023.112607).
+  [arXiv: 2211.14009](https://doi.org/10.48550/arXiv.2211.14009).
+  
+This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl.
+"""
+@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr,
+                                                                         normal_direction_avg::AbstractVector,
+                                                                         equations::IdealGlmMhdEquations3D,
+                                                                         nonconservative_type::NonConservativeSymmetric,
+                                                                         nonconservative_term::Integer)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
+
+    if nonconservative_term == 1
+        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+        # The factor 0.5 of the average can be omitted since it is already applied when this
+        # function is called.
+        B_dot_n_avg = ((B1_ll + B1_rr) * normal_direction_avg[1] +
+                       (B2_ll + B2_rr) * normal_direction_avg[2] +
+                       (B3_ll + B3_rr) * normal_direction_avg[3])
+        f = SVector(0,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    0)
+    else # nonconservative_term == 2
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        # The factor 0.5 of the average can be omitted since it is already applied when this
+        # function is called.
+        psi_avg = (psi_ll + psi_rr)
+        f = SVector(0,
+                    0,
+                    0,
+                    0,
+                    psi_avg,
+                    0,
+                    0,
+                    0,
+                    psi_avg)
+    end
+
+    return f
+end
+
 """
     flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations3D)
 
@@ -1130,6 +1313,20 @@ p &= (\gamma - 1) \left( E_\mathrm{tot} - E_\mathrm{kin} - E_\mathrm{mag} - \fra
     return p
 end
 
+# Transformation from conservative variables u to d(p)/d(u)
+@inline function gradient_conservative(::typeof(pressure),
+                                       u, equations::IdealGlmMhdEquations3D)
+    rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
+
+    v1 = rho_v1 / rho
+    v2 = rho_v2 / rho
+    v3 = rho_v3 / rho
+    v_square = v1^2 + v2^2 + v3^2
+
+    return (equations.gamma - 1) *
+           SVector(0.5f0 * v_square, -v1, -v2, -v3, 1, -B1, -B2, -B3, -psi)
+end
+
 @inline function density_pressure(u, equations::IdealGlmMhdEquations3D)
     rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3, psi = u
     rho_times_p = (equations.gamma - 1) * (rho * rho_e -
@@ -1407,6 +1604,14 @@ end
             cons[9]^2 / 2)
 end
 
+# State validation for Newton-bisection method of subcell IDP limiting
+@inline function Base.isvalid(u, equations::IdealGlmMhdEquations3D)
+    if u[1] <= 0 || pressure(u, equations) <= 0
+        return false
+    end
+    return true
+end
+
 # Calculate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons'
 @inline function cross_helicity(cons, ::IdealGlmMhdEquations3D)
     return (cons[2] * cons[6] + cons[3] * cons[7] + cons[4] * cons[8]) / cons[1]
diff --git a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl
index 7b10c3117a5..b7ef4c24c88 100644
--- a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl
+++ b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl
@@ -29,7 +29,26 @@ function create_cache(mesh::P4estMesh{3},
                                                                  nnodes(dg))
 
     if have_nonconservative_terms(equations) == true
-        error("Unsupported system of equations with nonconservative terms")
+        A5d = Array{uEltype, 5}
+        # Extract the nonconservative flux as a dispatch argument for `n_nonconservative_terms`
+        _, volume_flux_noncons = volume_integral.volume_flux_dg
+
+        flux_nonconservative_temp_threaded = A5d[A5d(undef, nvariables(equations),
+                                                     n_nonconservative_terms(volume_flux_noncons),
+                                                     nnodes(dg), nnodes(dg),
+                                                     nnodes(dg))
+                                                 for _ in 1:Threads.maxthreadid()]
+        fhat_nonconservative_temp_threaded = A5d[A5d(undef, nvariables(equations),
+                                                     n_nonconservative_terms(volume_flux_noncons),
+                                                     nnodes(dg), nnodes(dg),
+                                                     nnodes(dg))
+                                                 for _ in 1:Threads.maxthreadid()]
+        phi_threaded = A5d[A5d(undef, nvariables(equations),
+                               n_nonconservative_terms(volume_flux_noncons),
+                               nnodes(dg), nnodes(dg), nnodes(dg))
+                           for _ in 1:Threads.maxthreadid()]
+        cache = (; cache..., flux_nonconservative_temp_threaded,
+                 fhat_nonconservative_temp_threaded, phi_threaded)
     end
 
     return (; cache..., antidiffusive_fluxes,
@@ -231,6 +250,298 @@ end
     return nothing
 end
 
+# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element
+# (**with non-conservative terms in "local * symmetric" form**).
+#
+# See also `flux_differencing_kernel!`.
+#
+# The calculation of the non-conservative staggered "fluxes" requires non-conservative
+# terms that can be written as a product of local and a symmetric contributions. See, e.g.,
+#
+# - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+#   Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.
+#
+@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R,
+                                u, mesh::P4estMesh{3},
+                                nonconservative_terms::True, equations,
+                                volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM,
+                                element,
+                                cache) where {
+                                              F_CONS <: Function,
+                                              F_NONCONS <:
+                                              FluxNonConservative{NonConservativeSymmetric()}
+                                              }
+    (; contravariant_vectors) = cache.elements
+    (; weights, derivative_split) = dg.basis
+    (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache
+    (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache
+
+    volume_flux_cons, volume_flux_noncons = volume_flux
+
+    flux_temp = flux_temp_threaded[Threads.threadid()]
+    flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()]
+
+    fhat_temp = fhat_temp_threaded[Threads.threadid()]
+    fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()]
+    phi = phi_threaded[Threads.threadid()]
+
+    # The FV-form fluxes are calculated in a recursive manner, i.e.:
+    # fhat_(0,1)   = w_0 * FVol_0,
+    # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j,   for j=1,...,N-1,
+    # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i).
+
+    # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated
+    # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing`
+    # and saved in in `flux_temp`.
+
+    # Split form volume flux in orientation 1: x direction
+    flux_temp .= zero(eltype(flux_temp))
+    flux_noncons_temp .= zero(eltype(flux_noncons_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja1_node = get_contravariant_vector(1, contravariant_vectors,
+                                            i, j, k, element) # x direction
+
+        # All diagonal entries of `derivative_split` are zero. Thus, we can skip
+        # the computation of the diagonal terms. In addition, we use the symmetry
+        # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux
+        # computations.
+        for ii in (i + 1):nnodes(dg)
+            u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
+            # pull the contravariant vectors and compute the average
+            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors,
+                                                   ii, j, k, element)
+            Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
+
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde1 = volume_flux_cons(u_node, u_node_ii, Ja1_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,
+                                       equations, dg, ii, j, k)
+            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+                # We multiply by 0.5 because that is done in other parts of Trixi
+                flux1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_avg,
+                                                    equations,
+                                                    NonConservativeSymmetric(), noncons)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[i, ii],
+                                           flux1_noncons,
+                                           equations, dg, noncons, i, j, k)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[ii, i],
+                                           flux1_noncons,
+                                           equations, dg, noncons, ii, j, k)
+            end
+        end
+    end
+
+    # FV-form flux `fhat` in x direction
+    fhat_temp[:, 1, :, :] .= zero(eltype(fhat1_L))
+    fhat_noncons_temp[:, :, 1, :, :] .= zero(eltype(fhat1_L))
+
+    # Compute local contribution to non-conservative flux
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_local = get_node_vars(u, equations, dg, i, j, k, element)
+        # pull the local contravariant vector
+        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+            set_node_vars!(phi,
+                           volume_flux_noncons(u_local, Ja1_node, equations,
+                                               NonConservativeLocal(), noncons),
+                           equations, dg, noncons, i, j, k)
+        end
+    end
+
+    for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1)
+        # Conservative part
+        for v in eachvariable(equations)
+            value = fhat_temp[v, i, j, k] + weights[i] * flux_temp[v, i, j, k]
+            fhat_temp[v, i + 1, j, k] = value
+            fhat1_L[v, i + 1, j, k] = value
+            fhat1_R[v, i + 1, j, k] = value
+        end
+        # Nonconservative part
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),
+            v in eachvariable(equations)
+
+            value = fhat_noncons_temp[v, noncons, i, j, k] +
+                    weights[i] * flux_noncons_temp[v, noncons, i, j, k]
+            fhat_noncons_temp[v, noncons, i + 1, j, k] = value
+
+            fhat1_L[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] +
+                                      phi[v, noncons, i, j, k] * value
+            fhat1_R[v, i + 1, j, k] = fhat1_R[v, i + 1, j, k] +
+                                      phi[v, noncons, i + 1, j, k] * value
+        end
+    end
+
+    # Split form volume flux in orientation 2: y direction
+    flux_temp .= zero(eltype(flux_temp))
+    flux_noncons_temp .= zero(eltype(flux_noncons_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja2_node = get_contravariant_vector(2, contravariant_vectors,
+                                            i, j, k, element)
+
+        for jj in (j + 1):nnodes(dg)
+            u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
+            # pull the contravariant vectors and compute the average
+            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors,
+                                                   i, jj, k, element)
+            Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde2 = volume_flux_cons(u_node, u_node_jj, Ja2_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,
+                                       equations, dg, i, jj, k)
+            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+                # We multiply by 0.5 because that is done in other parts of Trixi
+                flux2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_avg,
+                                                    equations,
+                                                    NonConservativeSymmetric(), noncons)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[j, jj],
+                                           flux2_noncons,
+                                           equations, dg, noncons, i, j, k)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[jj, j],
+                                           flux2_noncons,
+                                           equations, dg, noncons, i, jj, k)
+            end
+        end
+    end
+
+    # FV-form flux `fhat` in y direction
+    fhat_temp[:, :, 1, :] .= zero(eltype(fhat1_L))
+    fhat_noncons_temp[:, :, :, 1, :] .= zero(eltype(fhat1_L))
+
+    # Compute local contribution to non-conservative flux
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_local = get_node_vars(u, equations, dg, i, j, k, element)
+        # pull the local contravariant vector
+        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+            set_node_vars!(phi,
+                           volume_flux_noncons(u_local, Ja2_node, equations,
+                                               NonConservativeLocal(), noncons),
+                           equations, dg, noncons, i, j, k)
+        end
+    end
+
+    for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg)
+        # Conservative part
+        for v in eachvariable(equations)
+            value = fhat_temp[v, i, j, k] + weights[j] * flux_temp[v, i, j, k]
+            fhat_temp[v, i, j + 1, k] = value
+            fhat2_L[v, i, j + 1, k] = value
+            fhat2_R[v, i, j + 1, k] = value
+        end
+        # Nonconservative part
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),
+            v in eachvariable(equations)
+
+            value = fhat_noncons_temp[v, noncons, i, j, k] +
+                    weights[j] * flux_noncons_temp[v, noncons, i, j, k]
+            fhat_noncons_temp[v, noncons, i, j + 1, k] = value
+
+            fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] +
+                                      phi[v, noncons, i, j, k] * value
+            fhat2_R[v, i, j + 1, k] = fhat2_R[v, i, j + 1, k] +
+                                      phi[v, noncons, i, j + 1, k] * value
+        end
+    end
+
+    # Split form volume flux in orientation 3: z direction
+    flux_temp .= zero(eltype(flux_temp))
+    flux_noncons_temp .= zero(eltype(flux_noncons_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja3_node = get_contravariant_vector(3, contravariant_vectors,
+                                            i, j, k, element)
+
+        for kk in (k + 1):nnodes(dg)
+            u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
+            # pull the contravariant vectors and compute the average
+            Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors,
+                                                   i, j, kk, element)
+            Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk)
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde3 = volume_flux_cons(u_node, u_node_kk, Ja3_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[k, kk], fluxtilde3,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[kk, k], fluxtilde3,
+                                       equations, dg, i, j, kk)
+            for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+                # We multiply by 0.5 because that is done in other parts of Trixi
+                flux3_noncons = volume_flux_noncons(u_node, u_node_kk, Ja3_avg,
+                                                    equations,
+                                                    NonConservativeSymmetric(), noncons)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[k, kk],
+                                           flux3_noncons,
+                                           equations, dg, noncons, i, j, k)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[kk, k],
+                                           flux3_noncons,
+                                           equations, dg, noncons, i, j, kk)
+            end
+        end
+    end
+
+    # FV-form flux `fhat` in z direction
+    fhat_temp[:, :, :, 1] .= zero(eltype(fhat1_L))
+    fhat_noncons_temp[:, :, :, :, 1] .= zero(eltype(fhat1_L))
+
+    # Compute local contribution to non-conservative flux
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_local = get_node_vars(u, equations, dg, i, j, k, element)
+        # pull the local contravariant vector
+        Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons)
+            set_node_vars!(phi,
+                           volume_flux_noncons(u_local, Ja3_node, equations,
+                                               NonConservativeLocal(), noncons),
+                           equations, dg, noncons, i, j, k)
+        end
+    end
+
+    for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg)
+        # Conservative part
+        for v in eachvariable(equations)
+            value = fhat_temp[v, i, j, k] + weights[k] * flux_temp[v, i, j, k]
+            fhat_temp[v, i, j, k + 1] = value
+            fhat3_L[v, i, j, k + 1] = value
+            fhat3_R[v, i, j, k + 1] = value
+        end
+        # Nonconservative part
+        for noncons in 1:n_nonconservative_terms(volume_flux_noncons),
+            v in eachvariable(equations)
+
+            value = fhat_noncons_temp[v, noncons, i, j, k] +
+                    weights[k] * flux_noncons_temp[v, noncons, i, j, k]
+            fhat_noncons_temp[v, noncons, i, j, k + 1] = value
+
+            fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] +
+                                      phi[v, noncons, i, j, k] * value
+            fhat3_R[v, i, j, k + 1] = fhat3_R[v, i, j, k + 1] +
+                                      phi[v, noncons, i, j, k + 1] * value
+        end
+    end
+
+    return nothing
+end
+
 # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems.
 @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R,
                                          fhat2_L, fhat2_R,
@@ -281,4 +592,52 @@ end
 
     return nothing
 end
+
+# Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems.
+@inline function calcflux_antidiffusive!(fhat1_L, fhat1_R,
+                                         fhat2_L, fhat2_R,
+                                         fhat3_L, fhat3_R,
+                                         fstar1_L, fstar1_R,
+                                         fstar2_L, fstar2_R,
+                                         fstar3_L, fstar3_R,
+                                         u, mesh::P4estMesh{3},
+                                         nonconservative_terms::True, equations,
+                                         limiter::SubcellLimiterIDP, dg, element, cache)
+    @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes
+
+    # Due to the use of LGL nodes, the DG staggered fluxes `fhat` and FV fluxes `fstar` are equal
+    # on the element interfaces. So, they are not computed in the volume integral and set to zero
+    # in their respective computation.
+    # The antidiffusive fluxes are therefore zero on the element interfaces and don't need to be
+    # computed either. They are set to zero directly after resizing the container.
+    # This applies to the indices `i=1` and `i=nnodes(dg)+1` for `antidiffusive_flux1_L` and
+    # `antidiffusive_flux1_R` and analogously for the other two directions.
+
+    for k in eachnode(dg), j in eachnode(dg), i in 2:nnodes(dg)
+        for v in eachvariable(equations)
+            antidiffusive_flux1_L[v, i, j, k, element] = fhat1_L[v, i, j, k] -
+                                                         fstar1_L[v, i, j, k]
+            antidiffusive_flux1_R[v, i, j, k, element] = fhat1_R[v, i, j, k] -
+                                                         fstar1_R[v, i, j, k]
+        end
+    end
+    for k in eachnode(dg), j in 2:nnodes(dg), i in eachnode(dg)
+        for v in eachvariable(equations)
+            antidiffusive_flux2_L[v, i, j, k, element] = fhat2_L[v, i, j, k] -
+                                                         fstar2_L[v, i, j, k]
+            antidiffusive_flux2_R[v, i, j, k, element] = fhat2_R[v, i, j, k] -
+                                                         fstar2_R[v, i, j, k]
+        end
+    end
+    for k in 2:nnodes(dg), j in eachnode(dg), i in eachnode(dg)
+        for v in eachvariable(equations)
+            antidiffusive_flux3_L[v, i, j, k, element] = fhat3_L[v, i, j, k] -
+                                                         fstar3_L[v, i, j, k]
+            antidiffusive_flux3_R[v, i, j, k, element] = fhat3_R[v, i, j, k] -
+                                                         fstar3_R[v, i, j, k]
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
index aed469a962a..0e63d83d8dd 100644
--- a/test/test_p4est_3d.jl
+++ b/test/test_p4est_3d.jl
@@ -661,6 +661,40 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_mhd_shockcapturing_subcell.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_subcell.jl"),
+                        l2=[
+                            0.006729931970167595,
+                            0.008638393158639436,
+                            0.008978257148101689,
+                            0.0085466685190268,
+                            0.0285664608641833,
+                            0.005796806835751598,
+                            0.007485378539184046,
+                            0.005846216235895686,
+                            1.116115482859488e-5
+                        ],
+                        linf=[
+                            0.31507940417652525,
+                            0.27581230560179737,
+                            0.5096527712168957,
+                            0.2900021150087706,
+                            0.9484970527977867,
+                            0.2591599747174065,
+                            0.22934145164154485,
+                            0.2868673643088755,
+                            0.0013663401454538622
+                        ],
+                        tspan=(0.0, 0.04))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    # Larger values for allowed allocations due to usage of custom
+    # integrator which are not *recorded* for the methods from
+    # OrdinaryDiffEq.jl
+    # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
+    @test_allocations(Trixi.rhs!, semi, sol, 15_000)
+end
+
 @trixi_testset "elixir_mhd_amr_entropy_bounded.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_amr_entropy_bounded.jl"),
                         l2=[