diff --git a/Project.toml b/Project.toml
index 296883f5420..95e1da07fd0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -108,7 +108,7 @@ StaticArrays = "1.5"
 StrideArrays = "0.1.26"
 StructArrays = "0.6.11"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.7.2"
+T8code = "0.7.4"
 TimerOutputs = "0.5.7"
 Triangulate = "2.2"
 TriplotBase = "0.1"
diff --git a/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl b/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl
new file mode 100644
index 00000000000..7efffbe52e3
--- /dev/null
+++ b/examples/t8code_3d_dgsem/elixir_advection_cubed_sphere.jl
@@ -0,0 +1,61 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+initial_condition = initial_condition_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = Dict(:inside => boundary_condition,
+                           :outside => boundary_condition)
+
+# Note that the first argument refers to the level of refinement, unlike in for p4est
+mesh = Trixi.T8codeMeshCubedSphere(5, 3, 0.5, 0.5;
+                                   polydeg = 3, initial_refinement_level = 0)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 100,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.2)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl
new file mode 100644
index 00000000000..d3509069681
--- /dev/null
+++ b/examples/t8code_3d_dgsem/elixir_euler_baroclinic_instability.jl
@@ -0,0 +1,299 @@
+# An idealized baroclinic instability test case
+# For optimal results consider increasing the resolution to 16x16x8 trees per cube face.
+#
+# Note that this elixir can take several hours to run.
+# Using 24 threads of an AMD Ryzen Threadripper 3990X (more threads don't speed it up further)
+# and `check-bounds=no`, this elixirs takes about one hour to run.
+# With 16x16x8 trees per cube face on the same machine, it takes about 28 hours.
+#
+# References:
+# - Paul A. Ullrich, Thomas Melvin, Christiane Jablonowski, Andrew Staniforth (2013)
+#   A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores
+#   https://doi.org/10.1002/qj.2241
+
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra
+
+###############################################################################
+# Setup for the baroclinic instability test
+gamma = 1.4
+equations = CompressibleEulerEquations3D(gamma)
+
+# Initial condition for an idealized baroclinic instability test
+# https://doi.org/10.1002/qj.2241, Section 3.2 and Appendix A
+function initial_condition_baroclinic_instability(x, t,
+                                                  equations::CompressibleEulerEquations3D)
+    lon, lat, r = cartesian_to_sphere(x)
+    radius_earth = 6.371229e6
+    # Make sure that the r is not smaller than radius_earth
+    z = max(r - radius_earth, 0.0)
+
+    # Unperturbed basic state
+    rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z)
+
+    # Stream function type perturbation
+    u_perturbation, v_perturbation = perturbation_stream_function(lon, lat, z)
+
+    u += u_perturbation
+    v = v_perturbation
+
+    # Convert spherical velocity to Cartesian
+    v1 = -sin(lon) * u - sin(lat) * cos(lon) * v
+    v2 = cos(lon) * u - sin(lat) * sin(lon) * v
+    v3 = cos(lat) * v
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+# Steady state for RHS correction below
+function steady_state_baroclinic_instability(x, t, equations::CompressibleEulerEquations3D)
+    lon, lat, r = cartesian_to_sphere(x)
+    radius_earth = 6.371229e6
+    # Make sure that the r is not smaller than radius_earth
+    z = max(r - radius_earth, 0.0)
+
+    # Unperturbed basic state
+    rho, u, p = basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z)
+
+    # Convert spherical velocity to Cartesian
+    v1 = -sin(lon) * u
+    v2 = cos(lon) * u
+    v3 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+function cartesian_to_sphere(x)
+    r = norm(x)
+    lambda = atan(x[2], x[1])
+    if lambda < 0
+        lambda += 2 * pi
+    end
+    phi = asin(x[3] / r)
+
+    return lambda, phi, r
+end
+
+# Unperturbed balanced steady-state.
+# Returns primitive variables with only the velocity in longitudinal direction (rho, u, p).
+# The other velocity components are zero.
+function basic_state_baroclinic_instability_longitudinal_velocity(lon, lat, z)
+    # Parameters from Table 1 in the paper
+    # Corresponding names in the paper are commented
+    radius_earth = 6.371229e6  # a
+    half_width_parameter = 2           # b
+    gravitational_acceleration = 9.80616     # g
+    k = 3           # k
+    surface_pressure = 1e5         # p₀
+    gas_constant = 287         # R
+    surface_equatorial_temperature = 310.0       # T₀ᴱ
+    surface_polar_temperature = 240.0       # T₀ᴾ
+    lapse_rate = 0.005       # Γ
+    angular_velocity = 7.29212e-5  # Ω
+
+    # Distance to the center of the Earth
+    r = z + radius_earth
+
+    # In the paper: T₀
+    temperature0 = 0.5 * (surface_equatorial_temperature + surface_polar_temperature)
+    # In the paper: A, B, C, H
+    const_a = 1 / lapse_rate
+    const_b = (temperature0 - surface_polar_temperature) /
+              (temperature0 * surface_polar_temperature)
+    const_c = 0.5 * (k + 2) * (surface_equatorial_temperature - surface_polar_temperature) /
+              (surface_equatorial_temperature * surface_polar_temperature)
+    const_h = gas_constant * temperature0 / gravitational_acceleration
+
+    # In the paper: (r - a) / bH
+    scaled_z = z / (half_width_parameter * const_h)
+
+    # Temporary variables
+    temp1 = exp(lapse_rate / temperature0 * z)
+    temp2 = exp(-scaled_z^2)
+
+    # In the paper: ̃τ₁, ̃τ₂
+    tau1 = const_a * lapse_rate / temperature0 * temp1 +
+           const_b * (1 - 2 * scaled_z^2) * temp2
+    tau2 = const_c * (1 - 2 * scaled_z^2) * temp2
+
+    # In the paper: ∫τ₁(r') dr', ∫τ₂(r') dr'
+    inttau1 = const_a * (temp1 - 1) + const_b * z * temp2
+    inttau2 = const_c * z * temp2
+
+    # Temporary variables
+    temp3 = r / radius_earth * cos(lat)
+    temp4 = temp3^k - k / (k + 2) * temp3^(k + 2)
+
+    # In the paper: T
+    temperature = 1 / ((r / radius_earth)^2 * (tau1 - tau2 * temp4))
+
+    # In the paper: U, u (zonal wind, first component of spherical velocity)
+    big_u = gravitational_acceleration / radius_earth * k * temperature * inttau2 *
+            (temp3^(k - 1) - temp3^(k + 1))
+    temp5 = radius_earth * cos(lat)
+    u = -angular_velocity * temp5 + sqrt(angular_velocity^2 * temp5^2 + temp5 * big_u)
+
+    # Hydrostatic pressure
+    p = surface_pressure *
+        exp(-gravitational_acceleration / gas_constant * (inttau1 - inttau2 * temp4))
+
+    # Density (via ideal gas law)
+    rho = p / (gas_constant * temperature)
+
+    return rho, u, p
+end
+
+# Perturbation as in Equations 25 and 26 of the paper (analytical derivative)
+function perturbation_stream_function(lon, lat, z)
+    # Parameters from Table 1 in the paper
+    # Corresponding names in the paper are commented
+    perturbation_radius = 1 / 6      # d₀ / a
+    perturbed_wind_amplitude = 1.0      # Vₚ
+    perturbation_lon = pi / 9     # Longitude of perturbation location
+    perturbation_lat = 2 * pi / 9 # Latitude of perturbation location
+    pertz = 15000    # Perturbation height cap
+
+    # Great circle distance (d in the paper) divided by a (radius of the Earth)
+    # because we never actually need d without dividing by a
+    great_circle_distance_by_a = acos(sin(perturbation_lat) * sin(lat) +
+                                      cos(perturbation_lat) * cos(lat) *
+                                      cos(lon - perturbation_lon))
+
+    # In the first case, the vertical taper function is per definition zero.
+    # In the second case, the stream function is per definition zero.
+    if z > pertz || great_circle_distance_by_a > perturbation_radius
+        return 0.0, 0.0
+    end
+
+    # Vertical tapering of stream function
+    perttaper = 1.0 - 3 * z^2 / pertz^2 + 2 * z^3 / pertz^3
+
+    # sin/cos(pi * d / (2 * d_0)) in the paper
+    sin_, cos_ = sincos(0.5 * pi * great_circle_distance_by_a / perturbation_radius)
+
+    # Common factor for both u and v
+    factor = 16 / (3 * sqrt(3)) * perturbed_wind_amplitude * perttaper * cos_^3 * sin_
+
+    u_perturbation = -factor * (-sin(perturbation_lat) * cos(lat) +
+                      cos(perturbation_lat) * sin(lat) * cos(lon - perturbation_lon)) /
+                     sin(great_circle_distance_by_a)
+
+    v_perturbation = factor * cos(perturbation_lat) * sin(lon - perturbation_lon) /
+                     sin(great_circle_distance_by_a)
+
+    return u_perturbation, v_perturbation
+end
+
+@inline function source_terms_baroclinic_instability(u, x, t,
+                                                     equations::CompressibleEulerEquations3D)
+    radius_earth = 6.371229e6  # a
+    gravitational_acceleration = 9.80616     # g
+    angular_velocity = 7.29212e-5  # Ω
+
+    r = norm(x)
+    # Make sure that r is not smaller than radius_earth
+    z = max(r - radius_earth, 0.0)
+    r = z + radius_earth
+
+    du1 = zero(eltype(u))
+
+    # Gravity term
+    temp = -gravitational_acceleration * radius_earth^2 / r^3
+    du2 = temp * u[1] * x[1]
+    du3 = temp * u[1] * x[2]
+    du4 = temp * u[1] * x[3]
+    du5 = temp * (u[2] * x[1] + u[3] * x[2] + u[4] * x[3])
+
+    # Coriolis term, -2Ω × ρv = -2 * angular_velocity * (0, 0, 1) × u[2:4]
+    du2 -= -2 * angular_velocity * u[3]
+    du3 -= 2 * angular_velocity * u[2]
+
+    return SVector(du1, du2, du3, du4, du5)
+end
+
+###############################################################################
+# Start of the actual elixir, semidiscretization of the problem
+
+initial_condition = initial_condition_baroclinic_instability
+
+boundary_conditions = Dict(:inside => boundary_condition_slip_wall,
+                           :outside => boundary_condition_slip_wall)
+
+# This is a good estimate for the speed of sound in this example.
+# Other values between 300 and 400 should work as well.
+surface_flux = FluxLMARS(340)
+volume_flux = flux_kennedy_gruber
+solver = DGSEM(polydeg = 5, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+# For optimal results, use (16, 8) here
+trees_per_cube_face = (8, 4)
+mesh = Trixi.T8codeMeshCubedSphere(trees_per_cube_face..., 6.371229e6, 30000.0,
+                                   polydeg = 5, initial_refinement_level = 0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_baroclinic_instability,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10 * 24 * 60 * 60.0) # time in seconds for 10 days
+
+# Save RHS of the steady state and subtract it in every RHS evaluation.
+# This trick preserves the steady state exactly (to machine rounding errors, of course).
+# Otherwise, this elixir produces entirely unusable results for a resolution of 8x8x4 cells
+# per cube face with a polydeg of 3.
+# With this trick, even the polydeg 3 simulation produces usable (although badly resolved) results,
+# and most of the grid imprinting in higher polydeg simulation is eliminated.
+#
+# See https://github.com/trixi-framework/Trixi.jl/issues/980 for more information.
+u_steady_state = compute_coefficients(steady_state_baroclinic_instability, tspan[1], semi)
+# Use a `let` block for performance (otherwise du_steady_state will be a global variable)
+let du_steady_state = similar(u_steady_state)
+    # Save RHS of the steady state
+    Trixi.rhs!(du_steady_state, u_steady_state, semi, tspan[1])
+
+    global function corrected_rhs!(du, u, semi, t)
+        # Normal RHS evaluation
+        Trixi.rhs!(du, u, semi, t)
+        # Correct by subtracting the steady-state RHS
+        Trixi.@trixi_timeit Trixi.timer() "rhs correction" begin
+            # Use Trixi.@threaded for threaded performance
+            Trixi.@threaded for i in eachindex(du)
+                du[i] -= du_steady_state[i]
+            end
+        end
+    end
+end
+u0 = compute_coefficients(tspan[1], semi)
+ode = ODEProblem(corrected_rhs!, u0, tspan, semi)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 5000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 5000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+
+# Use a Runge-Kutta method with automatic (error based) time step size control
+# Enable threading of the RK method for better performance on multiple threads
+sol = solve(ode, RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()); abstol = 1.0e-6,
+            reltol = 1.0e-6,
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary
diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl
index fa5b1d8f81a..4c85c543e16 100644
--- a/src/meshes/t8code_mesh.jl
+++ b/src/meshes/t8code_mesh.jl
@@ -598,7 +598,7 @@ struct AbaqusFile{NDIMS}
 end
 
 """
-    T8codeMesh(meshfile::String, NDIMS; kwargs...)
+    T8codeMesh(filepath::String, NDIMS; kwargs...)
 
 Main mesh constructor for the `T8codeMesh` that imports an unstructured, conforming
 mesh from either a Gmsh mesh file (`.msh`) or Abaqus mesh file  (`.inp`) which is determined
@@ -783,6 +783,52 @@ function t8_cmesh_new_from_connectivity(connectivity::Ptr{p8est_connectivity}, c
     return t8_cmesh_new_from_p8est(connectivity, comm, 0)
 end
 
+"""
+T8codeMeshCubedSphere(trees_per_face_dimension, layers, inner_radius, thickness;
+                      polydeg, RealT=Float64, initial_refinement_level=0)
+
+Construct a cubed spherical shell of given inner radius and thickness as `T8codeMesh` with
+`6 * trees_per_face_dimension^2 * layers` trees. The mesh will have two boundaries,
+`:inside` and `:outside`.
+
+# Arguments
+- `lat_lon_levels_per_face_dimension::Integer`: number of trees per patch in longitudinal
+                                                and latitudinal direction given as level of
+                                                refinement.
+- `layers::Integer`: the number of trees in the third local dimension of each face, i.e.,
+                     the number of layers of the shell.
+- `inner_radius::Float64`: Radius of the inner side of the shell.
+- `thickness::Float64`: Thickness of the shell. The outer radius will be
+                        `inner_radius + thickness`.
+- `polydeg::Integer`: polynomial degree used to store the geometry of the mesh.
+                      The mapping will be approximated by an interpolation polynomial
+                      of the specified degree for each tree.
+- `RealT::Type`: the type that should be used for coordinates.
+- `initial_refinement_level::Integer`: refine the mesh uniformly to this level before the
+                                       simulation starts.
+"""
+function T8codeMeshCubedSphere(lat_lon_levels_per_face_dimension, layers, inner_radius,
+                               thickness;
+                               polydeg, RealT = Float64, initial_refinement_level = 0)
+    NDIMS = 3
+    cmesh = t8_cmesh_new_cubed_spherical_shell(inner_radius, thickness,
+                                               lat_lon_levels_per_face_dimension,
+                                               layers, mpi_comm())
+    do_face_ghost = mpi_isparallel()
+    scheme = t8_scheme_new_default_cxx()
+    forest = t8_forest_new_uniform(cmesh, scheme, initial_refinement_level, do_face_ghost,
+                                   mpi_comm())
+
+    num_trees = t8_cmesh_get_num_trees(cmesh)
+    boundary_names = fill(Symbol("---"), 2 * NDIMS, num_trees)
+    for itree in 1:num_trees
+        boundary_names[5, itree] = :inside
+        boundary_names[6, itree] = :outside
+    end
+
+    return T8codeMesh{NDIMS, RealT}(forest, boundary_names; polydeg = polydeg)
+end
+
 struct adapt_callback_passthrough
     adapt_callback::Function
     user_data::Any
diff --git a/src/solvers/dgsem_t8code/containers_3d.jl b/src/solvers/dgsem_t8code/containers_3d.jl
index 981dabc397f..a4eb86ebad3 100644
--- a/src/solvers/dgsem_t8code/containers_3d.jl
+++ b/src/solvers/dgsem_t8code/containers_3d.jl
@@ -64,8 +64,7 @@ function calc_node_coordinates!(node_coordinates,
                                          current_index += 1),
                                     matrix1, matrix2, matrix3,
                                     view(mesh.tree_node_coordinates, :, :, :, :,
-                                         global_itree + 1),
-                                    tmp1)
+                                         global_itree + 1), tmp1)
         end
     end
 
diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl
index e07cba1f576..82bec769f28 100644
--- a/test/test_t8code_3d.jl
+++ b/test/test_t8code_3d.jl
@@ -111,6 +111,22 @@ mkdir(outdir)
         end
     end
 
+    # This test differs from the one in `test_p4est_3d.jl` in the latitudinal and
+    # longitudinal dimensions.
+    @trixi_testset "elixir_advection_cubed_sphere.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"),
+                            l2=[0.002006918015656413],
+                            linf=[0.027655117058380085])
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
+
     # This test is identical to the one in `test_p4est_3d.jl`.
     @trixi_testset "elixir_advection_restart.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
@@ -336,6 +352,38 @@ mkdir(outdir)
         end
     end
 
+    # This test is identical to the one in `test_p4est_3d.jl`.
+    @trixi_testset "elixir_euler_baroclinic_instability.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "elixir_euler_baroclinic_instability.jl"),
+                            l2=[
+                                6.725093801700048e-7,
+                                0.00021710076010951073,
+                                0.0004386796338203878,
+                                0.00020836270267103122,
+                                0.07601887903440395
+                            ],
+                            linf=[
+                                1.9107530539574924e-5,
+                                0.02980358831035801,
+                                0.048476331898047564,
+                                0.02200137344113612,
+                                4.848310144356219
+                            ],
+                            tspan=(0.0, 1e2),
+                            # Decrease tolerance of adaptive time stepping to get similar results across different systems
+                            abstol=1.0e-9, reltol=1.0e-9,
+                            coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+    end
+
     @trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
                             l2=[