diff --git a/examples/tree_2d_dgsem/elixir_euler_density_wave_nonconforming_idp_mortars.jl b/examples/tree_2d_dgsem/elixir_euler_density_wave_nonconforming_idp_mortars.jl
new file mode 100644
index 00000000000..2547c7486a5
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_density_wave_nonconforming_idp_mortars.jl
@@ -0,0 +1,76 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_density_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 5
+basis = LobattoLegendreBasis(polydeg)
+
+# No limiting is needed in the volume integral
+# It is still required to use `VolumeIntegralSubcellLimiting` in order to use `MortarIDP`.
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = [],
+                                positivity_variables_nonlinear = [])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+mortar = MortarIDP(equations, basis;
+                   positivity_variables_cons = ["rho"])
+solver = DGSEM(basis, surface_flux, volume_integral, mortar)
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+refinement_patches = ((type = "box", coordinates_min = (-0.5, -0.5),
+                       coordinates_max = (0.0, 1.0)),
+                      (type = "box", coordinates_min = (0.5, -1.0),
+                       coordinates_max = (1.0, 0.5)))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                refinement_patches = refinement_patches,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim,
+                                     extra_node_variables = (:limiting_coefficient,))
+
+stepsize_callback = StepsizeCallback(cfl = 0.4)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_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/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl
similarity index 74%
rename from examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl
rename to examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl
index 5d47dec8f96..5bbd268f28b 100644
--- a/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl
@@ -30,25 +30,20 @@ function initial_condition_kelvin_helmholtz_instability(x, t,
 end
 initial_condition = initial_condition_kelvin_helmholtz_instability
 
-# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
-# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
-# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
-# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
-# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
-# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the 
-# `StepsizeCallback` (CFL-Condition) and less diffusion.
-surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
+surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 positivity_variables_cons = ["rho"],
-                                positivity_variables_nonlinear = [pressure])
+                                positivity_variables_nonlinear = [pressure],
+                                local_twosided_variables_cons = []) # required for testing
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+mortar = MortarIDP(equations, basis; positivity_variables_cons = ["rho"])
+solver = DGSEM(basis, surface_flux, volume_integral, mortar)
 
 coordinates_min = (-1.0, -1.0)
 coordinates_max = (1.0, 1.0)
@@ -79,12 +74,27 @@ save_solution = SaveSolutionCallback(interval = 100,
 save_restart = SaveRestartCallback(interval = 1000,
                                    save_final_restart = true)
 
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max = 1.0,
+                                          alpha_min = 0.0001,
+                                          alpha_smooth = false,
+                                          variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 4,
+                                      med_level = 0, med_threshold = 0.0003,
+                                      max_level = 6, max_threshold = 0.003)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 1,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        stepsize_callback,
-                        save_restart, save_solution)
+                        amr_callback,
+                        save_restart, save_solution,
+                        stepsize_callback)
 
 ###############################################################################
 # run the simulation
@@ -95,5 +105,4 @@ 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
-                  ode_default_options()...,
                   callback = callbacks);
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 6b596bd74fb..964de6f25e8 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -265,7 +265,7 @@ export DG,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
-       MortarL2
+       MortarL2, MortarIDP
 
 export reconstruction_O2_inner, reconstruction_O2_full,
        reconstruction_constant,
diff --git a/src/callbacks_stage/subcell_limiter_idp_correction.jl b/src/callbacks_stage/subcell_limiter_idp_correction.jl
index e1cb42035d1..ebe5f6747dd 100644
--- a/src/callbacks_stage/subcell_limiter_idp_correction.jl
+++ b/src/callbacks_stage/subcell_limiter_idp_correction.jl
@@ -56,6 +56,19 @@ function (limiter!::SubcellLimiterIDPCorrection)(u_ode, semi, t, dt,
 
     perform_idp_correction!(u, dt, mesh, equations, solver, cache)
 
+    if solver.mortar isa Trixi.LobattoLegendreMortarIDP
+        @trixi_timeit timer() "mortar blending factors" calc_mortar_limiting_factor!(u,
+                                                                                     semi,
+                                                                                     t,
+                                                                                     dt)
+
+        @trixi_timeit timer() "mortar correction" perform_idp_mortar_correction(u, dt,
+                                                                                mesh,
+                                                                                equations,
+                                                                                solver,
+                                                                                cache)
+    end
+
     return nothing
 end
 
diff --git a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl
index 4caaff8fc17..b93a62c5d8d 100644
--- a/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl
+++ b/src/callbacks_stage/subcell_limiter_idp_correction_2d.jl
@@ -67,4 +67,112 @@ function perform_idp_correction!(u, dt,
 
     return nothing
 end
+
+function perform_idp_mortar_correction(u, dt, mesh::TreeMesh{2}, equations, dg, cache)
+    (; orientations, limiting_factor) = cache.mortars
+
+    (; surface_flux_values) = cache.elements
+    (; surface_flux_values_high_order) = cache.antidiffusive_fluxes
+    (; boundary_interpolation) = dg.basis
+
+    for mortar in eachmortar(dg, cache)
+        if isapprox(limiting_factor[mortar], one(eltype(limiting_factor)))
+            continue
+        end
+
+        if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+            if orientations[mortar] == 1
+                direction_small = 1
+                direction_large = 2
+            else
+                direction_small = 3
+                direction_large = 4
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_small = boundary_interpolation[1, 1]
+            factor_large = -boundary_interpolation[nnodes(dg), 2]
+        else # large_sides[mortar] == 2 -> small elements on left side
+            if orientations[mortar] == 1
+                direction_small = 2
+                direction_large = 1
+            else
+                direction_small = 4
+                direction_large = 3
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_large = boundary_interpolation[1, 1]
+            factor_small = -boundary_interpolation[nnodes(dg), 2]
+        end
+
+        for i in eachnode(dg)
+            if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (1, i)
+                    indices_large = (nnodes(dg), i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, 1)
+                    indices_large = (i, nnodes(dg))
+                end
+            else # large_sides[mortar] == 2 -> small elements on left side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (nnodes(dg), i)
+                    indices_large = (1, i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, nnodes(dg))
+                    indices_large = (i, 1)
+                end
+            end
+
+            # small elements
+            for small_element_index in 1:2
+                small_element = cache.mortars.neighbor_ids[small_element_index, mortar]
+                inverse_jacobian_small = get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                              mesh, indices_small...,
+                                                              small_element)
+
+                flux_small_high_order = get_node_vars(surface_flux_values_high_order,
+                                                      equations, dg,
+                                                      i, direction_small, small_element)
+                flux_small_low_order = get_node_vars(surface_flux_values, equations, dg,
+                                                     i, direction_small, small_element)
+                flux_difference_small = factor_small *
+                                        (flux_small_high_order .- flux_small_low_order)
+
+                multiply_add_to_node_vars!(u,
+                                           dt * inverse_jacobian_small *
+                                           (1 - limiting_factor[mortar]),
+                                           flux_difference_small, equations, dg,
+                                           indices_small..., small_element)
+            end
+
+            # large element
+            large_element = cache.mortars.neighbor_ids[3, mortar]
+            inverse_jacobian_large = get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                          mesh, indices_large...,
+                                                          large_element)
+
+            flux_large_high_order = get_node_vars(surface_flux_values_high_order,
+                                                  equations, dg, i, direction_large,
+                                                  large_element)
+            flux_large_low_order = get_node_vars(surface_flux_values, equations, dg, i,
+                                                 direction_large, large_element)
+            flux_difference_large = factor_large *
+                                    (flux_large_high_order .- flux_large_low_order)
+
+            multiply_add_to_node_vars!(u,
+                                       dt * inverse_jacobian_large *
+                                       (1 - limiting_factor[mortar]),
+                                       flux_difference_large, equations, dg,
+                                       indices_large..., large_element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd
diff --git a/src/solvers/dgsem/basis_lobatto_legendre.jl b/src/solvers/dgsem/basis_lobatto_legendre.jl
index 9647f172e20..a62b94f7f72 100644
--- a/src/solvers/dgsem/basis_lobatto_legendre.jl
+++ b/src/solvers/dgsem/basis_lobatto_legendre.jl
@@ -77,13 +77,13 @@ function LobattoLegendreBasis(RealT, polydeg::Integer)
     derivative_split_transpose = Matrix(derivative_split')
     derivative_dhat = calc_dhat(nodes_, weights_)
 
-    # Type conversions to enable possible optimizations of runtime performance 
+    # Type conversions to enable possible optimizations of runtime performance
     # and latency
     nodes = SVector{nnodes_, RealT}(nodes_)
     weights = SVector{nnodes_, RealT}(weights_)
     inverse_weights = SVector{nnodes_, RealT}(inverse_weights_)
 
-    # We keep the matrices above stored using the standard `Matrix` type 
+    # We keep the matrices above stored using the standard `Matrix` type
     # since this is usually as fast as `SMatrix`
     # (when using `let` in the volume integral/`@threaded`)
     # and reduces latency
@@ -139,7 +139,7 @@ end
     eachnode(basis::LobattoLegendreBasis)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the nodes in `basis`. 
+for the nodes in `basis`.
 In particular, not the nodes themselves are returned.
 """
 @inline eachnode(basis::LobattoLegendreBasis) = Base.OneTo(nnodes(basis))
@@ -201,7 +201,7 @@ function MortarL2(basis::LobattoLegendreBasis)
     reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
     reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
 
-    # We keep the matrices above stored using the standard `Matrix` type 
+    # We keep the matrices above stored using the standard `Matrix` type
     # since this is usually as fast as `SMatrix`
     # (when using `let` in the volume integral/`@threaded`)
     # and reduces latency
@@ -242,6 +242,60 @@ end
 
 @inline polydeg(mortar::LobattoLegendreMortarL2) = nnodes(mortar) - 1
 
+struct LobattoLegendreMortarIDP{RealT <: Real, NNODES, NDIMS, LENGTH, Mortar} <:
+       AbstractMortar{RealT}
+    positivity_variables_cons::Vector{Int}
+    mortar_l2::Mortar
+    # LENGTH = `2 * (NDIMS - 1) + 1`
+    mortar_weights::Array{RealT, LENGTH}     # [node_i (large), node_j (large), node_i (small), node_j (small), small_element]
+    mortar_weights_sums::Array{RealT, NDIMS} # [node_i, node_j, small (1) / large (2) element]
+end
+
+function MortarIDP(equations, basis::LobattoLegendreBasis;
+                   positivity_variables_cons = String[])
+    RealT = real(basis)
+    n_dims = ndims(equations)
+    nnodes_ = nnodes(basis)
+
+    mortar_l2 = MortarL2(basis)
+
+    mortar_weights, mortar_weights_sums = calc_mortar_weights(equations, basis, RealT)
+
+    positivity_variables_cons_ = get_variable_index.(positivity_variables_cons,
+                                                     equations)
+
+    LobattoLegendreMortarIDP{RealT, nnodes_, n_dims,
+                             2 * (n_dims - 1) + 1, typeof(mortar_l2)}(positivity_variables_cons_,
+                                                                      mortar_l2,
+                                                                      mortar_weights,
+                                                                      mortar_weights_sums)
+end
+
+function Base.show(io::IO, mortar::LobattoLegendreMortarIDP)
+    @nospecialize mortar # reduce precompilation time
+
+    print(io, "LobattoLegendreMortarIDP(polydeg=", polydeg(mortar), ", positivity for ",
+          mortar.positivity_variables_cons, ")")
+end
+function Base.show(io::IO, ::MIME"text/plain", mortar::LobattoLegendreMortarIDP)
+    @nospecialize mortar # reduce precompilation time
+
+    print(io, "LobattoLegendreMortarIDP{", real(mortar),
+          "} with polynomials of degree ",
+          polydeg(mortar),
+          " and positivity limiting for conservative variables ",
+          mortar.positivity_variables_cons)
+end
+
+@inline Base.real(mortar::LobattoLegendreMortarIDP{RealT}) where {RealT} = RealT
+
+@inline function nnodes(mortar::LobattoLegendreMortarIDP{RealT, NNODES}) where {RealT,
+                                                                                NNODES}
+    NNODES
+end
+
+@inline polydeg(mortar::LobattoLegendreMortarIDP) = nnodes(mortar) - 1
+
 # TODO: We can create EC mortars along the lines of the following implementation.
 # abstract type AbstractMortarEC{RealT} <: AbstractMortar{RealT} end
 
@@ -292,7 +346,7 @@ function SolutionAnalyzer(basis::LobattoLegendreBasis;
     nodes_, weights_ = gauss_lobatto_nodes_weights(nnodes_, RealT)
     vandermonde = polynomial_interpolation_matrix(get_nodes(basis), nodes_)
 
-    # Type conversions to enable possible optimizations of runtime performance 
+    # Type conversions to enable possible optimizations of runtime performance
     # and latency
     nodes = SVector{nnodes_, RealT}(nodes_)
     weights = SVector{nnodes_, RealT}(weights_)
@@ -325,7 +379,7 @@ end
     eachnode(analyzer::LobattoLegendreAnalyzer)
 
 Return an iterator over the indices that specify the location in relevant data structures
-for the nodes in `analyzer`. 
+for the nodes in `analyzer`.
 In particular, not the nodes themselves are returned.
 """
 @inline eachnode(analyzer::LobattoLegendreAnalyzer) = Base.OneTo(nnodes(analyzer))
@@ -531,19 +585,19 @@ function calc_lhat(x, nodes, weights)
     return lhat
 end
 
-""" 
+"""
     lagrange_interpolating_polynomials(x, nodes, wbary)
 
 Calculate Lagrange polynomials for a given node distribution with
-associated barycentric weights `wbary` at a given point `x` on the 
+associated barycentric weights `wbary` at a given point `x` on the
 reference interval ``[-1, 1]``.
 
 This returns all ``l_j(x)``, i.e., the Lagrange polynomials for each node ``x_j``.
-Thus, to obtain the interpolating polynomial ``p(x)`` at ``x``, one has to 
+Thus, to obtain the interpolating polynomial ``p(x)`` at ``x``, one has to
 multiply the Lagrange polynomials with the nodal values ``u_j`` and sum them up:
 ``p(x) = \\sum_{j=1}^{n} u_j l_j(x)``.
 
-For details, see e.g. Section 2 of 
+For details, see e.g. Section 2 of
 - Jean-Paul Berrut and Lloyd N. Trefethen (2004).
   Barycentric Lagrange Interpolation.
   [DOI:10.1137/S0036144502417715](https://doi.org/10.1137/S0036144502417715)
@@ -553,7 +607,7 @@ function lagrange_interpolating_polynomials(x, nodes, wbary)
     polynomials = zeros(eltype(nodes), n_nodes)
 
     for i in 1:n_nodes
-        # Avoid division by zero when `x` is close to node by using 
+        # Avoid division by zero when `x` is close to node by using
         # the Kronecker-delta property at nodes
         # of the Lagrange interpolation polynomials.
         if isapprox(x, nodes[i], rtol = eps(x))
@@ -580,9 +634,9 @@ end
 Computes nodes ``x_j`` and weights ``w_j`` for the (Legendre-)Gauss-Lobatto quadrature.
 This implements algorithm 25 "GaussLobattoNodesAndWeights" from the book
 
-- David A. Kopriva, (2009). 
+- David A. Kopriva, (2009).
   Implementing spectral methods for partial differential equations:
-  Algorithms for scientists and engineers. 
+  Algorithms for scientists and engineers.
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
 function gauss_lobatto_nodes_weights(n_nodes::Integer, RealT = Float64)
@@ -686,9 +740,9 @@ end
 Computes nodes ``x_j`` and weights ``w_j`` for the Gauss-Legendre quadrature.
 This implements algorithm 23 "LegendreGaussNodesAndWeights" from the book
 
-- David A. Kopriva, (2009). 
+- David A. Kopriva, (2009).
   Implementing spectral methods for partial differential equations:
-  Algorithms for scientists and engineers. 
+  Algorithms for scientists and engineers.
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
 function gauss_nodes_weights(n_nodes::Integer, RealT = Float64)
@@ -756,9 +810,9 @@ end
 Computes the Legendre polynomial of degree `N` and its derivative at `x`.
 This implements algorithm 22 "LegendrePolynomialAndDerivative" from the book
 
-- David A. Kopriva, (2009). 
+- David A. Kopriva, (2009).
   Implementing spectral methods for partial differential equations:
-  Algorithms for scientists and engineers. 
+  Algorithms for scientists and engineers.
   [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
 """
 function legendre_polynomial_and_derivative(N::Int, x::Real)
diff --git a/src/solvers/dgsem_tree/containers_2d.jl b/src/solvers/dgsem_tree/containers_2d.jl
index f7467c2db48..4c480ee337a 100644
--- a/src/solvers/dgsem_tree/containers_2d.jl
+++ b/src/solvers/dgsem_tree/containers_2d.jl
@@ -605,6 +605,139 @@ function init_mortars(cell_ids, mesh::TreeMesh2D,
     return mortars
 end
 
+# Container data structure (structure-of-arrays style) for DG mortars for IDP AMR
+# Positions/directions for orientations = 1, large_sides = 2:
+# mortar is orthogonal to x-axis, large side is in positive coordinate direction wrt mortar
+#           |    |
+# upper = 2 |    |
+#           |    |
+#                | 3 = large side
+#           |    |
+# lower = 1 |    |
+#           |    |
+mutable struct IDPMortarContainer2D{uEltype <: Real} <: AbstractContainer
+    u_upper::Array{uEltype, 4}  # [leftright, variables, i, mortars]
+    u_lower::Array{uEltype, 4}  # [leftright, variables, i, mortars]
+    u_large::Array{uEltype, 3}  # [variables, i, mortars]
+    neighbor_ids::Array{Int, 2} # [position, mortars]
+    # Large sides: left -> 1, right -> 2
+    large_sides::Vector{Int}  # [mortars]
+    orientations::Vector{Int} # [mortars]
+    limiting_factor::Vector{uEltype} # [mortars]
+    # internal `resize!`able storage
+    _u_upper::Vector{uEltype}
+    _u_lower::Vector{uEltype}
+    _u_large::Vector{uEltype}
+    _neighbor_ids::Vector{Int}
+end
+
+nvariables(mortars::IDPMortarContainer2D) = size(mortars.u_upper, 2)
+nnodes(mortars::IDPMortarContainer2D) = size(mortars.u_upper, 3)
+Base.eltype(mortars::IDPMortarContainer2D) = eltype(mortars.u_upper)
+
+# See explanation of Base.resize! for the element container
+function Base.resize!(mortars::IDPMortarContainer2D, capacity)
+    n_nodes = nnodes(mortars)
+    n_variables = nvariables(mortars)
+    @unpack _u_upper, _u_lower, _u_large, _neighbor_ids,
+    large_sides, orientations, limiting_factor = mortars
+
+    resize!(_u_upper, 2 * n_variables * n_nodes * capacity)
+    mortars.u_upper = unsafe_wrap(Array, pointer(_u_upper),
+                                  (2, n_variables, n_nodes, capacity))
+
+    resize!(_u_lower, 2 * n_variables * n_nodes * capacity)
+    mortars.u_lower = unsafe_wrap(Array, pointer(_u_lower),
+                                  (2, n_variables, n_nodes, capacity))
+
+    resize!(_u_large, n_variables * n_nodes * capacity)
+    mortars.u_large = unsafe_wrap(Array, pointer(_u_large),
+                                  (n_variables, n_nodes, capacity))
+
+    resize!(_neighbor_ids, 3 * capacity)
+    mortars.neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids),
+                                       (3, capacity))
+
+    resize!(large_sides, capacity)
+
+    resize!(orientations, capacity)
+
+    resize!(limiting_factor, capacity)
+
+    return nothing
+end
+
+function IDPMortarContainer2D{uEltype}(capacity::Integer, n_variables,
+                                       n_nodes) where {uEltype <: Real}
+    nan = convert(uEltype, NaN)
+
+    # Initialize fields with defaults
+    _u_upper = fill(nan, 2 * n_variables * n_nodes * capacity)
+    u_upper = unsafe_wrap(Array, pointer(_u_upper),
+                          (2, n_variables, n_nodes, capacity))
+
+    _u_lower = fill(nan, 2 * n_variables * n_nodes * capacity)
+    u_lower = unsafe_wrap(Array, pointer(_u_lower),
+                          (2, n_variables, n_nodes, capacity))
+
+    _u_large = fill(nan, n_variables * n_nodes * capacity)
+    u_large = unsafe_wrap(Array, pointer(_u_large),
+                          (n_variables, n_nodes, capacity))
+
+    _neighbor_ids = fill(typemin(Int), 3 * capacity)
+    neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids),
+                               (3, capacity))
+
+    large_sides = fill(typemin(Int), capacity)
+
+    orientations = fill(typemin(Int), capacity)
+
+    limiting_factor = fill(nan, capacity)
+
+    return IDPMortarContainer2D{uEltype}(u_upper, u_lower, u_large, neighbor_ids,
+                                         large_sides, orientations, limiting_factor,
+                                         _u_upper, _u_lower, _u_large, _neighbor_ids)
+end
+
+# Return number of IDP mortars
+@inline nmortars(l2mortars::IDPMortarContainer2D) = length(l2mortars.orientations)
+
+# Allow printing container contents
+function Base.show(io::IO, ::MIME"text/plain", c::IDPMortarContainer2D)
+    @nospecialize c # reduce precompilation time
+
+    println(io, '*'^20)
+    for idx in CartesianIndices(c.u_upper)
+        println(io, "c.u_upper[$idx] = $(c.u_upper[idx])")
+    end
+    for idx in CartesianIndices(c.u_lower)
+        println(io, "c.u_lower[$idx] = $(c.u_lower[idx])")
+    end
+    for idx in CartesianIndices(c.u_large)
+        println(io, "c.u_large[$idx] = $(c.u_large[idx])")
+    end
+    println(io, "transpose(c.neighbor_ids) = $(transpose(c.neighbor_ids))")
+    println(io, "c.large_sides = $(c.large_sides)")
+    println(io, "c.orientations = $(c.orientations)")
+    println(io, "c.limiting_factor = $(c.limiting_factor)")
+    print(io, '*'^20)
+end
+
+# Create mortar container and initialize mortar data in `elements`.
+function init_mortars(cell_ids, mesh::TreeMesh2D,
+                      elements::ElementContainer2D,
+                      mortar::LobattoLegendreMortarIDP)
+    # Initialize containers
+    n_mortars = count_required_mortars(mesh, cell_ids)
+    mortars = IDPMortarContainer2D{eltype(elements)}(n_mortars,
+                                                     nvariables(elements),
+                                                     nnodes(elements))
+
+    # Connect elements with mortars
+    init_mortars!(mortars, elements, mesh)
+    return mortars
+end
+
 # Count the number of mortars that need to be created
 function count_required_mortars(mesh::TreeMesh2D, cell_ids)
     count = 0
@@ -1272,11 +1405,13 @@ mutable struct ContainerAntidiffusiveFlux2D{uEltype <: Real} <: AbstractContaine
     antidiffusive_flux1_R::Array{uEltype, 4} # [variables, i, j, elements]
     antidiffusive_flux2_L::Array{uEltype, 4} # [variables, i, j, elements]
     antidiffusive_flux2_R::Array{uEltype, 4} # [variables, i, j, elements]
+    surface_flux_values_high_order::Array{uEltype, 4} # [variables, i, direction, elements]
     # internal `resize!`able storage
     _antidiffusive_flux1_L::Vector{uEltype}
     _antidiffusive_flux1_R::Vector{uEltype}
     _antidiffusive_flux2_L::Vector{uEltype}
     _antidiffusive_flux2_R::Vector{uEltype}
+    _surface_flux_values_high_order::Vector{uEltype}
 end
 
 function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables,
@@ -1302,14 +1437,23 @@ function ContainerAntidiffusiveFlux2D{uEltype}(capacity::Integer, n_variables,
     antidiffusive_flux2_R = unsafe_wrap(Array, pointer(_antidiffusive_flux2_R),
                                         (n_variables, n_nodes, n_nodes + 1, capacity))
 
+    _surface_flux_values_high_order = fill(nan_uEltype,
+                                           n_variables * n_nodes * 2 * 2 * capacity)
+    surface_flux_values_high_order = unsafe_wrap(Array,
+                                                 pointer(_surface_flux_values_high_order),
+                                                 (n_variables, n_nodes, 2 * 2,
+                                                  capacity))
+
     return ContainerAntidiffusiveFlux2D{uEltype}(antidiffusive_flux1_L,
                                                  antidiffusive_flux1_R,
                                                  antidiffusive_flux2_L,
                                                  antidiffusive_flux2_R,
+                                                 surface_flux_values_high_order,
                                                  _antidiffusive_flux1_L,
                                                  _antidiffusive_flux1_R,
                                                  _antidiffusive_flux2_L,
-                                                 _antidiffusive_flux2_R)
+                                                 _antidiffusive_flux2_R,
+                                                 _surface_flux_values_high_order)
 end
 
 nvariables(fluxes::ContainerAntidiffusiveFlux2D) = size(fluxes.antidiffusive_flux1_L, 1)
@@ -1324,7 +1468,7 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity)
     n_nodes = nnodes(fluxes)
     n_variables = nvariables(fluxes)
 
-    @unpack _antidiffusive_flux1_L, _antidiffusive_flux2_L, _antidiffusive_flux1_R, _antidiffusive_flux2_R = fluxes
+    @unpack _antidiffusive_flux1_L, _antidiffusive_flux2_L, _antidiffusive_flux1_R, _antidiffusive_flux2_R, _surface_flux_values_high_order = fluxes
 
     resize!(_antidiffusive_flux1_L, n_variables * (n_nodes + 1) * n_nodes * capacity)
     fluxes.antidiffusive_flux1_L = unsafe_wrap(Array, pointer(_antidiffusive_flux1_L),
@@ -1343,6 +1487,12 @@ function Base.resize!(fluxes::ContainerAntidiffusiveFlux2D, capacity)
                                                (n_variables, n_nodes, n_nodes + 1,
                                                 capacity))
 
+    resize!(_surface_flux_values_high_order, n_variables * n_nodes * 2 * 2 * capacity)
+    fluxes.surface_flux_values_high_order = unsafe_wrap(Array,
+                                                        pointer(_surface_flux_values_high_order),
+                                                        (n_variables, n_nodes, 2 * 2,
+                                                         capacity))
+
     uEltype = eltype(fluxes.antidiffusive_flux1_L)
     @threaded for element in axes(fluxes.antidiffusive_flux1_L, 4)
         fluxes.antidiffusive_flux1_L[:, 1, :, element] .= zero(uEltype)
diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 844d934294f..c424db3d707 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -86,7 +86,8 @@ end
 # The methods below are specialized on the mortar type
 # and called from the basic `create_cache` method at the top.
 function create_cache(mesh::TreeMesh{2}, equations,
-                      mortar_l2::LobattoLegendreMortarL2, uEltype)
+                      mortar_l2::Union{LobattoLegendreMortarL2,
+                                       LobattoLegendreMortarIDP}, uEltype)
     # TODO: Taal performance using different types
     MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)},
                   uEltype, 2,
diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl
index 11230746494..b4eccb55a82 100644
--- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl
+++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl
@@ -65,6 +65,73 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, P4estMesh{2}},
             flux_temp_threaded, fhat_temp_threaded)
 end
 
+function calc_mortar_weights(equations::AbstractEquations{2},
+                             basis::LobattoLegendreBasis, RealT)
+    n_nodes = nnodes(basis)
+    mortar_weights = zeros(RealT, n_nodes, n_nodes, 2) # [node_i (large element), node_i (small element), small element]
+    mortar_weights_sums = zeros(RealT, n_nodes, 2)     # [node, left (=1) or large (=2) element]
+
+    calc_mortar_weights!(equations, mortar_weights, n_nodes, RealT)
+
+    # Sums of mortar weights for normalization
+    for i in eachnode(basis)
+        for k in eachnode(basis)
+            # Add weights from large element to small element
+            # Sums for both small elements are equal due to symmetry
+            mortar_weights_sums[i, 1] += mortar_weights[k, i, 1]
+            # Add weights from small element to large element
+            for small_element in 1:2
+                mortar_weights_sums[i, 2] += mortar_weights[i, k, small_element]
+            end
+        end
+    end
+
+    return mortar_weights, mortar_weights_sums
+end
+
+function calc_mortar_weights!(equations::AbstractEquations{2},
+                              mortar_weights, n_nodes, RealT)
+    _, weights = gauss_lobatto_nodes_weights(n_nodes, RealT)
+
+    # Local mortar weights are of the form: `w_ij = int_S psi_i phi_j ds`,
+    # where `psi_i` are the basis functions of the large element and `phi_j` are the basis
+    # functions of the small element. `S` is the face connecting both elements.
+    # We use piecewise constant basis functions on the LGL subgrid. So, only focus on interval,
+    # where both basis functions are non-zero. `interval = [left_bound, right_bound]`.
+    # `w_ij = int_S psi_i phi_j ds = int_{left_bound}^{right_bound} ds = right_bound - left_bound`.
+    # `right_bound = min(left_bound_large, left_bound_small)`
+    # `left_bound = max(right_bound_large, right_bound_small)`
+    # If `right_bound <= left_bound`, i.e., both intervals don't overlap, then `w_ij = 0`.
+
+    # Due to the LGL subgrid, the interval bounds are cumulative LGL quadrature weights.
+    cum_weights_large = [zero(RealT); cumsum(weights)] .- 1 # on [-1, 1]
+    cum_weights_lower = 0.5f0 * cum_weights_large .- 0.5f0  # on [-1, 0]
+    cum_weights_upper = cum_weights_lower .+ 1              # on [0, 1]
+    # So, for `w_ij` we have
+    # `right_bound = min(cum_weights_large[i], cum_weights_small[j])`
+    # `left_bound = max(cum_weights_large[i+1], cum_weights_small[j+1])`
+
+    for j in 1:n_nodes, i in 1:n_nodes
+        # lower and large element element
+        left = max(cum_weights_large[i], cum_weights_lower[j])
+        right = min(cum_weights_large[i + 1], cum_weights_lower[j + 1])
+
+        # Local weight of 0 if intervals do not overlap, i.e., `right <= left`
+        if right > left
+            mortar_weights[i, j, 1] = right - left
+        end
+
+        # upper and large element
+        left = max(cum_weights_large[i], cum_weights_upper[j])
+        right = min(cum_weights_large[i + 1], cum_weights_upper[j + 1])
+        if right > left
+            mortar_weights[i, j, 2] = right - left
+        end
+    end
+
+    return mortar_weights
+end
+
 # Subcell limiting currently only implemented for certain mesh types
 function calc_volume_integral!(du, u,
                                mesh::Union{TreeMesh{2}, StructuredMesh{2},
@@ -722,6 +789,197 @@ end
     return nothing
 end
 
+function prolong2mortars!(cache, u, mesh::TreeMesh{2}, equations,
+                          mortar_idp::LobattoLegendreMortarIDP, dg::DGSEM)
+    prolong2mortars!(cache, u, mesh, equations, mortar_idp.mortar_l2, dg)
+
+    # The data of both small elements were already copied to the mortar cache
+    @threaded for mortar in eachmortar(dg, cache)
+        large_element = cache.mortars.neighbor_ids[3, mortar]
+
+        # Copy solutions
+        if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+            if cache.mortars.orientations[mortar] == 1
+                # IDP mortars in x-direction
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        cache.mortars.u_large[v, l, mortar] = u[v, nnodes(dg), l,
+                                                                large_element]
+                    end
+                end
+            else
+                # IDP mortars in y-direction
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        cache.mortars.u_large[v, l, mortar] = u[v, l, nnodes(dg),
+                                                                large_element]
+                    end
+                end
+            end
+        else # large_sides[mortar] == 2 -> small elements on left side
+            if cache.mortars.orientations[mortar] == 1
+                # IDP mortars in x-direction
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        cache.mortars.u_large[v, l, mortar] = u[v, 1, l, large_element]
+                    end
+                end
+            else
+                # IDP mortars in y-direction
+                for l in eachnode(dg)
+                    for v in eachvariable(equations)
+                        cache.mortars.u_large[v, l, mortar] = u[v, l, 1, large_element]
+                    end
+                end
+            end
+        end
+    end
+
+    return nothing
+end
+
+function calc_mortar_flux!(surface_flux_values, mesh,
+                           nonconservative_terms, equations,
+                           mortar_idp::LobattoLegendreMortarIDP, surface_integral,
+                           dg::DG, cache)
+    # low order fluxes
+    @trixi_timeit timer() "calc_mortar_flux_low_order!" calc_mortar_flux_low_order!(surface_flux_values,
+                                                                                    mesh,
+                                                                                    nonconservative_terms,
+                                                                                    equations,
+                                                                                    mortar_idp,
+                                                                                    surface_integral,
+                                                                                    dg,
+                                                                                    cache)
+
+    # high order fluxes
+    (; surface_flux_values_high_order) = cache.antidiffusive_fluxes
+    @trixi_timeit timer() "calc_mortar_flux!" calc_mortar_flux!(surface_flux_values_high_order,
+                                                                mesh,
+                                                                nonconservative_terms,
+                                                                equations,
+                                                                mortar_idp.mortar_l2,
+                                                                dg.surface_integral, dg,
+                                                                cache)
+
+    return nothing
+end
+
+function calc_mortar_flux_low_order!(surface_flux_values,
+                                     mesh::TreeMesh{2},
+                                     nonconservative_terms::False, equations,
+                                     mortar_idp::LobattoLegendreMortarIDP,
+                                     surface_integral, dg::DG, cache)
+    @unpack surface_flux = surface_integral
+    @unpack u_lower, u_upper, u_large, orientations = cache.mortars
+    (; mortar_weights, mortar_weights_sums) = mortar_idp
+
+    @threaded for mortar in eachmortar(dg, cache)
+        large_element = cache.mortars.neighbor_ids[3, mortar]
+        upper_element = cache.mortars.neighbor_ids[2, mortar]
+        lower_element = cache.mortars.neighbor_ids[1, mortar]
+
+        # Calculate fluxes
+        orientation = orientations[mortar]
+
+        if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+            if orientation == 1
+                # L2 mortars in x-direction
+                direction_small = 1
+                direction_large = 2
+            else
+                # L2 mortars in y-direction
+                direction_small = 3
+                direction_large = 4
+            end
+            small_side = 2
+        else # large_sides[mortar] == 2 -> small elements on left side
+            if orientation == 1
+                # L2 mortars in x-direction
+                direction_small = 2
+                direction_large = 1
+            else
+                # L2 mortars in y-direction
+                direction_small = 4
+                direction_large = 3
+            end
+            small_side = 1
+        end
+
+        surface_flux_values[:, :, direction_small, lower_element] .= zero(eltype(surface_flux_values))
+        surface_flux_values[:, :, direction_small, upper_element] .= zero(eltype(surface_flux_values))
+        surface_flux_values[:, :, direction_large, large_element] .= zero(eltype(surface_flux_values))
+        # Lower element
+        for i in eachnode(dg)
+            u_lower_local = get_surface_node_vars(u_lower, equations, dg,
+                                                  i, mortar)[small_side]
+            for k in eachnode(dg)
+                factor = mortar_weights[k, i, 1]
+                if isapprox(factor, zero(typeof(factor)))
+                    continue
+                end
+                u_large_local = get_node_vars(u_large, equations, dg, k, mortar)
+
+                if small_side == 2 # -> small elements on right side
+                    flux = surface_flux(u_large_local, u_lower_local, orientation,
+                                        equations)
+                else # small_side == 1 -> small elements on left side
+                    flux = surface_flux(u_lower_local, u_large_local, orientation,
+                                        equations)
+                end
+
+                # Lower element
+                multiply_add_to_node_vars!(surface_flux_values,
+                                           factor /
+                                           mortar_weights_sums[i, 1],
+                                           flux, equations, dg,
+                                           i, direction_small, lower_element)
+                # Large element
+                multiply_add_to_node_vars!(surface_flux_values,
+                                           factor /
+                                           mortar_weights_sums[k, 2],
+                                           flux, equations, dg,
+                                           k, direction_large, large_element)
+            end
+        end
+        # Upper element
+        for i in eachnode(dg)
+            u_upper_local = get_surface_node_vars(u_upper, equations, dg,
+                                                  i, mortar)[small_side]
+            for k in eachnode(dg)
+                factor = mortar_weights[k, i, 2]
+                if isapprox(factor, zero(typeof(factor)))
+                    continue
+                end
+                u_large_local = get_node_vars(u_large, equations, dg, k, mortar)
+
+                if small_side == 2 # -> small elements on right side
+                    flux = surface_flux(u_large_local, u_upper_local, orientation,
+                                        equations)
+                else # small_side == 1 -> small elements on left side
+                    flux = surface_flux(u_upper_local, u_large_local, orientation,
+                                        equations)
+                end
+
+                # Upper element
+                multiply_add_to_node_vars!(surface_flux_values,
+                                           factor /
+                                           mortar_weights_sums[i, 1],
+                                           flux, equations, dg,
+                                           i, direction_small, upper_element)
+                # Large element
+                multiply_add_to_node_vars!(surface_flux_values,
+                                           factor /
+                                           mortar_weights_sums[k, 2],
+                                           flux, equations, dg,
+                                           k, direction_large, large_element)
+            end
+        end
+    end
+
+    return nothing
+end
+
 """
     get_boundary_outer_state(u_inner, t,
                              boundary_condition::BoundaryConditionDirichlet,
diff --git a/src/solvers/dgsem_tree/subcell_limiters_2d.jl b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
index cca91aa94b0..ad93d4e7b73 100644
--- a/src/solvers/dgsem_tree/subcell_limiters_2d.jl
+++ b/src/solvers/dgsem_tree/subcell_limiters_2d.jl
@@ -83,6 +83,103 @@ end
         end
     end
 
+    # Calc bounds at mortars
+    # TODO: How to include values at mortar interfaces?
+    # - For LobattoLegendreMortarIDP: include only values of nodes with nonnegative local weights
+    # - For LobattoLegendreMortarL2: include all neighboring values (TODO?)
+    l2_mortars = dg.mortar isa LobattoLegendreMortarL2
+    for mortar in eachmortar(dg, cache)
+        large_element = cache.mortars.neighbor_ids[3, mortar]
+
+        orientation = cache.mortars.orientations[mortar]
+
+        for i in eachnode(dg)
+            if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+                if orientation == 1
+                    # L2 mortars in x-direction
+                    indices_small = (1, i)
+                    indices_large = (nnodes(dg), i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, 1)
+                    indices_large = (i, nnodes(dg))
+                end
+            else # large_sides[mortar] == 2 -> small elements on left side
+                if orientation == 1
+                    # L2 mortars in x-direction
+                    indices_small = (nnodes(dg), i)
+                    indices_large = (1, i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, nnodes(dg))
+                    indices_large = (i, 1)
+                end
+            end
+            # Get solution data
+            var_small = (u[variable, indices_small...,
+                           cache.mortars.neighbor_ids[1, mortar]],
+                         u[variable, indices_small...,
+                           cache.mortars.neighbor_ids[2, mortar]])
+            # Using the following version with `ntuple` creates allocations due to a type instability of `indices_small`.
+            # var_small = index -> u[variable, indices_small..., cache.mortars.neighbor_ids[index, mortar]]
+            # Theoretically, that could be fixed with the following version:
+            # f = let indices_small = indices_small
+            #     index -> u[variable, indices_small..., cache.mortars.neighbor_ids[index, mortar]]
+            # end
+            # var_small = ntuple(f, Val(2))
+            var_large = u[variable, indices_large..., large_element]
+
+            for j in eachnode(dg)
+                if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+                    if orientation == 1
+                        # L2 mortars in x-direction
+                        indices_small_inner = (1, j)
+                        indices_large_inner = (nnodes(dg), j)
+                    else
+                        # L2 mortars in y-direction
+                        indices_small_inner = (j, 1)
+                        indices_large_inner = (j, nnodes(dg))
+                    end
+                else # large_sides[mortar] == 2 -> small elements on left side
+                    if orientation == 1
+                        # L2 mortars in x-direction
+                        indices_small_inner = (nnodes(dg), j)
+                        indices_large_inner = (1, j)
+                    else
+                        # L2 mortars in y-direction
+                        indices_small_inner = (j, nnodes(dg))
+                        indices_large_inner = (j, 1)
+                    end
+                end
+
+                for small_element_index in 1:2
+                    small_element = cache.mortars.neighbor_ids[small_element_index,
+                                                               mortar]
+                    # from large to small element
+                    if l2_mortars ||
+                       dg.mortar.mortar_weights[i, j, small_element_index] > 0
+                        var_min[indices_small_inner..., small_element] = min(var_min[indices_small_inner...,
+                                                                                     small_element],
+                                                                             var_large)
+                        var_max[indices_small_inner..., small_element] = max(var_max[indices_small_inner...,
+                                                                                     small_element],
+                                                                             var_large)
+                    end
+                    # from small to large element
+                    if l2_mortars ||
+                       dg.mortar.mortar_weights[j, i, small_element_index] > 0
+                        var_min[indices_large_inner..., large_element] = min(var_min[indices_large_inner...,
+                                                                                     large_element],
+                                                                             var_small[small_element_index])
+                        var_max[indices_large_inner..., large_element] = max(var_max[indices_large_inner...,
+                                                                                     large_element],
+                                                                             var_small[small_element_index])
+                    end
+                end
+            end
+        end
+    end
+
     # Calc bounds at physical boundaries
     for boundary in eachboundary(dg, cache)
         element = cache.boundaries.neighbor_ids[boundary]
@@ -419,6 +516,214 @@ end
     return nothing
 end
 
+###############################################################################
+# IDP mortar limiting
+###############################################################################
+
+@inline function calc_mortar_limiting_factor!(u, semi, t, dt)
+    mesh, _, solver, cache = mesh_equations_solver_cache(semi)
+    (; positivity_variables_cons) = solver.mortar
+    (; limiting_factor) = cache.mortars
+    @trixi_timeit timer() "reset alpha" limiting_factor.=zeros(eltype(limiting_factor))
+
+    @trixi_timeit timer() "conservative variables" for var_index in positivity_variables_cons
+        limiting_positivity_conservative!(limiting_factor, u, dt, semi, mesh, var_index)
+    end
+
+    return nothing
+end
+
+###############################################################################
+# Local two-sided limiting of conservative variables
+@inline function limiting_positivity_conservative!(limiting_factor, u, dt, semi,
+                                                   mesh::TreeMesh{2}, var_index)
+    _, _, dg, cache = mesh_equations_solver_cache(semi)
+
+    (; orientations, large_sides, u_large, u_upper, u_lower) = cache.mortars
+    (; surface_flux_values, inverse_jacobian) = cache.elements
+    (; surface_flux_values_high_order) = cache.antidiffusive_fluxes
+    (; boundary_interpolation) = dg.basis
+
+    (; positivity_correction_factor) = dg.volume_integral.limiter
+
+    for mortar in eachmortar(dg, cache)
+        large_element = cache.mortars.neighbor_ids[3, mortar]
+        upper_element = cache.mortars.neighbor_ids[2, mortar]
+        lower_element = cache.mortars.neighbor_ids[1, mortar]
+
+        # Get small side of mortar
+        small_side = large_sides[mortar] == 1 ? 2 : 1
+
+        # Compute minimal bound
+        var_min_upper = typemax(eltype(surface_flux_values))
+        var_min_lower = typemax(eltype(surface_flux_values))
+        var_min_large = typemax(eltype(surface_flux_values))
+        for i in eachnode(dg)
+            if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (1, i)
+                    indices_large = (nnodes(dg), i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, 1)
+                    indices_large = (i, nnodes(dg))
+                end
+            else # large_sides[mortar] == 2 -> small elements on left side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (nnodes(dg), i)
+                    indices_large = (1, i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, nnodes(dg))
+                    indices_large = (i, 1)
+                end
+            end
+            var_upper = u[var_index, indices_small..., upper_element]
+            var_lower = u[var_index, indices_small..., lower_element]
+            var_large = u[var_index, indices_large..., large_element]
+            var_min_upper = min(var_min_upper, var_upper)
+            var_min_lower = min(var_min_lower, var_lower)
+            var_min_large = min(var_min_large, var_large)
+        end
+        var_min_upper = positivity_correction_factor * var_min_upper
+        var_min_lower = positivity_correction_factor * var_min_lower
+        var_min_large = positivity_correction_factor * var_min_large
+
+        # Set up correct direction and factors
+        if large_sides[mortar] == 1 # -> small elements on right side
+            if orientations[mortar] == 1
+                direction_small = 1
+                direction_large = 2
+            else
+                direction_small = 3
+                direction_large = 4
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_small = boundary_interpolation[1, 1]
+            factor_large = -boundary_interpolation[nnodes(dg), 2]
+        else # large_sides[mortar] == 2 -> small elements on left side
+            if orientations[mortar] == 1
+                direction_small = 2
+                direction_large = 1
+            else
+                direction_small = 4
+                direction_large = 3
+            end
+            # In `apply_jacobian`, `du` is multiplied with inverse jacobian and a negative sign.
+            # This sign switch is directly applied to the boundary interpolation factors here.
+            factor_large = boundary_interpolation[1, 1]
+            factor_small = -boundary_interpolation[nnodes(dg), 2]
+        end
+
+        # Compute limiting factor
+        for i in eachnode(dg)
+            if large_sides[mortar] == 1 # -> small elements on right side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (1, i)
+                    indices_large = (nnodes(dg), i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, 1)
+                    indices_large = (i, nnodes(dg))
+                end
+            else # large_sides[mortar] == 2 -> small elements on left side
+                if orientations[mortar] == 1
+                    # L2 mortars in x-direction
+                    indices_small = (nnodes(dg), i)
+                    indices_large = (1, i)
+                else
+                    # L2 mortars in y-direction
+                    indices_small = (i, nnodes(dg))
+                    indices_large = (i, 1)
+                end
+            end
+            var_upper = u[var_index, indices_small..., upper_element]
+            var_lower = u[var_index, indices_small..., lower_element]
+            var_large = u[var_index, indices_large..., large_element]
+
+            if min(var_upper, var_lower, var_large) < 0
+                error("Safe low-order method produces negative value for conservative variable rho. Try a smaller time step.")
+            end
+
+            # Compute flux differences
+            flux_lower_high_order = surface_flux_values_high_order[var_index, i,
+                                                                   direction_small,
+                                                                   lower_element]
+            flux_lower_low_order = surface_flux_values[var_index, i, direction_small,
+                                                       lower_element]
+            flux_difference_lower = factor_small *
+                                    (flux_lower_high_order - flux_lower_low_order)
+
+            flux_upper_high_order = surface_flux_values_high_order[var_index, i,
+                                                                   direction_small,
+                                                                   upper_element]
+            flux_upper_low_order = surface_flux_values[var_index, i, direction_small,
+                                                       upper_element]
+            flux_difference_upper = factor_small *
+                                    (flux_upper_high_order - flux_upper_low_order)
+
+            flux_large_high_order = surface_flux_values_high_order[var_index, i,
+                                                                   direction_large,
+                                                                   large_element]
+            flux_large_low_order = surface_flux_values[var_index, i, direction_large,
+                                                       large_element]
+            flux_difference_large = factor_large *
+                                    (flux_large_high_order - flux_large_low_order)
+
+            # Use pure low-order fluxes if high-order fluxes are not finite.
+            if !isfinite(flux_lower_high_order) ||
+               !isfinite(flux_upper_high_order) ||
+               !isfinite(flux_large_high_order)
+                limiting_factor[mortar] = 1
+                break
+            end
+
+            inverse_jacobian_upper = get_inverse_jacobian(inverse_jacobian, mesh,
+                                                          indices_small...,
+                                                          upper_element)
+            inverse_jacobian_lower = get_inverse_jacobian(inverse_jacobian, mesh,
+                                                          indices_small...,
+                                                          lower_element)
+            inverse_jacobian_large = get_inverse_jacobian(inverse_jacobian, mesh,
+                                                          indices_large...,
+                                                          large_element)
+
+            # Real one-sided Zalesak-type limiter
+            # * Zalesak (1979). "Fully multidimensional flux-corrected transport algorithms for fluids"
+            # * Kuzmin et al. (2010). "Failsafe flux limiting and constrained data projections for equations of gas dynamics"
+            # Note: The Zalesak limiter has to be computed, even if the state is valid, because the correction is
+            #       for each mortar, not each node
+            Qm_upper = min(0, var_min_upper - var_upper)
+            Qm_lower = min(0, var_min_lower - var_lower)
+            Qm_large = min(0, var_min_large - var_large)
+
+            Pm_upper = min(0, flux_difference_upper)
+            Pm_lower = min(0, flux_difference_lower)
+            Pm_large = min(0, flux_difference_large)
+
+            Pm_upper = dt * inverse_jacobian_upper * Pm_upper
+            Pm_lower = dt * inverse_jacobian_lower * Pm_lower
+            Pm_large = dt * inverse_jacobian_large * Pm_large
+
+            # Compute blending coefficient avoiding division by zero
+            # (as in paper of [Guermond, Nazarov, Popov, Thomas] (4.8))
+            Qm_upper = abs(Qm_upper) / (abs(Pm_upper) + eps(typeof(Qm_upper)) * 100)
+            Qm_lower = abs(Qm_lower) / (abs(Pm_lower) + eps(typeof(Qm_lower)) * 100)
+            Qm_large = abs(Qm_large) / (abs(Pm_large) + eps(typeof(Qm_large)) * 100)
+
+            # Calculate limiting factor
+            limiting_factor[mortar] = max(limiting_factor[mortar],
+                                          1 - Qm_upper, 1 - Qm_lower, 1 - Qm_large)
+        end
+    end
+
+    return nothing
+end
+
 ###############################################################################
 # Newton-bisection method
 
@@ -427,7 +732,8 @@ end
                                      variable, min_or_max,
                                      initial_check, final_check,
                                      inverse_jacobian, dt,
-                                     equations, dg, cache, limiter)
+                                     equations::AbstractEquations{2},
+                                     dg, cache, limiter)
     (; inverse_weights) = dg.basis
     (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
 
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index da5305f0217..8fc43209fa3 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -68,6 +68,32 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_euler_density_wave_nonconforming_idp_mortars.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_density_wave_nonconforming_idp_mortars.jl"),
+                        initial_refinement_level=2,
+                        l2=[
+                            0.01853467680355577,
+                            0.0018534676803570467,
+                            0.003706935360709665,
+                            0.0004633669200885599
+                        ],
+                        linf=[
+                            0.3234934074357616,
+                            0.03234934074339482,
+                            0.06469868148701247,
+                            0.008087335185543054
+                        ],
+                        tspan=(0.0, 0.5))
+    # 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_euler_source_terms_nonperiodic.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_source_terms_nonperiodic.jl"),
@@ -534,28 +560,30 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
-@trixi_testset "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl" begin
+@trixi_testset "elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl" begin
     rm(joinpath("out", "deviations.txt"), force = true)
     @test_trixi_include(joinpath(EXAMPLES_DIR,
-                                 "elixir_euler_kelvin_helmholtz_instability_sc_subcell.jl"),
+                                 "elixir_euler_kelvin_helmholtz_instability_amr_sc_subcell.jl"),
+                        local_twosided_variables_cons=["rho"],
+                        cfl=0.5,
                         l2=[
-                            0.42185634563805724,
-                            0.1686471269704017,
-                            0.18240674916968103,
-                            0.17858250604280654
+                            0.19657723128273924,
+                            0.09819990707800491,
+                            0.14969902881019737,
+                            0.07474283620648917
                         ],
                         linf=[
-                            1.7012978064377158,
-                            0.7149714986746726,
-                            0.5822547982757897,
-                            0.7300051017382696
+                            0.9087650823657859,
+                            0.48203285905738114,
+                            0.35330695561397685,
+                            0.3443844432700125
                         ],
-                        tspan=(0.0, 2.0),
+                        tspan=(0.0, 1.0),
                         save_errors=true)
     lines = readlines(joinpath("out", "deviations.txt"))
-    @test lines[1] == "# iter, simu_time, rho_min, pressure_min"
-    # Run takes 745 time steps
-    @test startswith(lines[end], "745")
+    @test lines[1] == "# iter, simu_time, rho_min, rho_max, pressure_min"
+    # Run takes 795 time steps
+    @test startswith(lines[end], "795")
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     # Larger values for allowed allocations due to usage of custom
@@ -563,6 +591,9 @@ end
     # OrdinaryDiffEq.jl
     # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
     @test_allocations(Trixi.rhs!, semi, sol, 15000)
+
+    # test long printing format
+    @test_nowarn display(solver.mortar)
 end
 
 @trixi_testset "elixir_euler_colliding_flow.jl" begin
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 037334fb617..daab704c329 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -452,6 +452,11 @@ end
     @test isnothing(display(c3d))
 end
 
+@timed_testset "DG IDP mortar container debug output" begin
+    c2d = Trixi.IDPMortarContainer2D{Float64}(1, 1, 1)
+    @test isnothing(display(c2d))
+end
+
 @timed_testset "TreeBoundaryContainer1D nnodes" begin
     capacity = 42
     n_variables = 9