diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 45a7f91..0e48049 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -37,11 +37,13 @@ jobs:
           # - 'nightly'
         os:
           - ubuntu-latest
-          - macOS-latest
           - windows-latest
         arch:
           - x64
         include:
+          - version: '1.9'
+            os: macOS-latest
+            arch: arm64
           - version: '1.7'
             os: ubuntu-latest
             arch: x64
diff --git a/Project.toml b/Project.toml
index f4a6db2..f8a6cd2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi2Vtk"
 uuid = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
 authors = ["Michael Schlottke-Lakemper <michael@sloede.com>", "Hendrik Ranocha <mail@ranocha.de>"]
-version = "0.3.16-pre"
+version = "0.3.19-DEV"
 
 [deps]
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
@@ -18,6 +18,6 @@ HDF5 = "0.14, 0.15, 0.16, 0.17"
 ProgressMeter = "1.3"
 StaticArrays = "0.12, 1.0"
 TimerOutputs = "0.5"
-Trixi = "0.5, 0.6, 0.7"
+Trixi = "0.5, 0.6, 0.7, 0.8, 0.9"
 WriteVTK = "1.7"
 julia = "1.7"
diff --git a/src/convert.jl b/src/convert.jl
index a66265d..c96dbe3 100644
--- a/src/convert.jl
+++ b/src/convert.jl
@@ -86,6 +86,10 @@ function trixi2vtk(filename::AbstractString...;
                         barlen = 40)
   end
 
+  # Show warning when reinterpolating node-level data of subcell limiting
+  # Auxiliary variable to show warning only once
+  has_warned_about_interpolation = false
+
   # Iterate over input files
   for (index, filename) in enumerate(filenames)
     verbose && println("Processing file $filename ($(index)/$(length(filenames)))...")
@@ -124,7 +128,7 @@ function trixi2vtk(filename::AbstractString...;
     if is_datafile
       verbose && println("| Reading data file...")
       @timeit "read data" (labels, data, n_elements, n_nodes,
-                           element_variables, time) = read_datafile(filename)
+                           element_variables, node_variables, time) = read_datafile(filename)
 
       assert_cells_elements(n_elements, mesh, filename, meshfile)
 
@@ -202,6 +206,28 @@ function trixi2vtk(filename::AbstractString...;
             verbose && println("| | Element variable: $label...")
             @timeit label vtk_celldata[label] = variable
           end
+
+          # Add node variables
+          for (label, variable) in node_variables
+            verbose && println("| | Node variable: $label...")
+            if reinterpolate
+              # Show warning if node-level data of subcell limiting are reinterpolated.
+              if label == "limiting_coefficient" && !has_warned_about_interpolation
+                println("WARNING: The limiting coefficients are no continuous field but happens " *
+                "to be represented by a piecewise-constant approximation. Thus, reinterpolation " *
+                "does not give a meaningful representation.")
+                has_warned_about_interpolation = true
+              end
+              @timeit "interpolate data" interpolated_cell_data = interpolate_data(Val(format),
+                                                                    reshape(variable, size(variable)..., 1),
+                                                                    mesh, n_visnodes, verbose)
+            else
+              @timeit "interpolate data" interpolated_cell_data = reshape(variable,
+                                                                          n_visnodes^ndims_ * n_elements)
+            end
+            # Add to node_data
+            @timeit label vtk_nodedata[label] = interpolated_cell_data
+          end
         end
       end
     end
diff --git a/src/io.jl b/src/io.jl
index 54b79fd..997cdf5 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -76,6 +76,21 @@ function read_datafile(filename::String)
       index +=1
     end
 
-    return labels, data, n_elements, n_nodes, element_variables, time
+    # Extract node variable arrays
+    node_variables = Dict{String, Union{Array{Float64}, Array{Int}}}()
+    index = 1
+    while haskey(file, "node_variables_$index")
+      varname = read(attributes(file["node_variables_$index"])["name"])
+      nodedata = read(file["node_variables_$index"])
+      if ndims_ == 2
+        node_variables[varname] = Array{Float64}(undef, n_nodes, n_nodes, n_elements)
+        @views node_variables[varname][:, :, :] .= nodedata
+      else
+        error("Unsupported number of spatial dimensions: ", ndims_)
+      end
+      index +=1
+    end
+
+    return labels, data, n_elements, n_nodes, element_variables, node_variables, time
   end
 end
diff --git a/test/Project.toml b/test/Project.toml
index 5d5e50e..adb8222 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,12 +2,15 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
+UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [compat]
 Documenter = "0.27, 1"
+Downloads = "1"
 OrdinaryDiffEq = "6"
 ReadVTK = "0.1, 0.2"
-Trixi = "0.5, 0.6, 0.7"
+Trixi = "0.5, 0.6, 0.7, 0.8, 0.9"
diff --git a/test/test_2d.jl b/test/test_2d.jl
index b661d0a..cd7ac31 100644
--- a/test/test_2d.jl
+++ b/test/test_2d.jl
@@ -22,8 +22,8 @@ end
 
     @timed_testset "mesh data" begin
       # create the output file to be tested
-      @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000010.h5"), output_directory=outdir)
-      outfilename = "mesh_000010_celldata.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "mesh_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir)
+      outfilename = "mesh_" * LEADING_ZEROS * "000010_celldata.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -38,8 +38,8 @@ end
 
     @timed_testset "solution celldata" begin
       # create the output file to be tested
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir)
-      outfilename = "solution_000010_celldata.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir)
+      outfilename = "solution_" * LEADING_ZEROS * "000010_celldata.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -54,8 +54,8 @@ end
 
     @timed_testset "reinterpolate with nonuniform data with VTU format" begin
       # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir)
-      outfilename = "solution_000010.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir)
+      outfilename = "solution_" * LEADING_ZEROS * "000010.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -70,8 +70,8 @@ end
 
     @timed_testset "reinterpolate with nonuniform data with VTI format" begin
       # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, format=:vti)
-      outfilename = "solution_000010.vti"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, format=:vti)
+      outfilename = "solution_" * LEADING_ZEROS * "000010.vti"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -86,8 +86,8 @@ end
 
     @timed_testset "do not reinterpolate with nonuniform data" begin
       # Create and test output without reinterpolation on LGL nodes
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false)
-      outfilename = "solution_000010.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=false)
+      outfilename = "solution_" * LEADING_ZEROS * "000010.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -103,8 +103,8 @@ end
     @timed_testset "do not reinterpolate with uniform data" begin
       # Create and test output without reinterpolation on uniform nodes
       # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-      outfilename = "solution_000010.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+      outfilename = "solution_" * LEADING_ZEROS * "000010.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -120,7 +120,7 @@ end
     @timed_testset "attempt reinterpolate with uniform data" begin
       # Purposely request a bad configuration and check that an error message gets thrown
       # OBS! Only needs tested once across all mesh types and dimensions
-      @test_throws ArgumentError trixi2vtk(joinpath(outdir, "solution_000010.h5"), output_directory=outdir, data_is_uniform=true)
+      @test_throws ArgumentError trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, data_is_uniform=true)
     end
   end
 
@@ -154,7 +154,7 @@ end
         #      that only needs to be tested once.
         @test_nowarn trixi2vtk(joinpath(outdir, "solution_00000*.h5"), output_directory=outdir)
 
-        outfilename = "solution_000001_celldata.vtu"
+        outfilename = "solution_" * LEADING_ZEROS * "000001_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -169,8 +169,8 @@ end
 
       @timed_testset "reinterpolate with nonuniform data" begin
         # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -185,8 +185,8 @@ end
 
       @timed_testset "do not reinterpolate with nonuniform data" begin
         # Create and test output without reinterpolation on LGL nodes
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -202,8 +202,8 @@ end
       @timed_testset "do not reinterpolate with uniform data" begin
         # Create and test output without reinterpolation on uniform nodes
         # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -219,7 +219,7 @@ end
       @timed_testset "attempt VTI format on unsupported mesh type" begin
         # Purposely request a bad configuration and check that an error message gets thrown
         # OBS! Only needs tested once across all mesh types and dimensions
-        @test_throws ArgumentError trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, format=:vti)
+        @test_throws ArgumentError trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, format=:vti)
       end
     end
 
@@ -245,8 +245,8 @@ end
 
       @timed_testset "solution celldata" begin
         # create the output file to be tested
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir)
-        outfilename = "solution_000001_celldata.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000001_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -261,8 +261,8 @@ end
 
       @timed_testset "reinterpolate with nonuniform data" begin
         # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -277,8 +277,8 @@ end
 
       @timed_testset "do not reinterpolate with nonuniform data" begin
         # Create and test output without reinterpolation on LGL nodes
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -294,8 +294,8 @@ end
       @timed_testset "do not reinterpolate with uniform data" begin
         # Create and test output without reinterpolation on uniform nodes
         # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -315,8 +315,8 @@ end
 
       @timed_testset "mesh data" begin
         # create the output file to be tested
-        @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000005.h5"), output_directory=outdir)
-        outfilename = "mesh_000005_celldata.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "mesh_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "mesh_" * LEADING_ZEROS * "000005_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -331,8 +331,8 @@ end
 
       @timed_testset "solution celldata" begin
         # create the output file to be tested
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir)
-        outfilename = "solution_000005_celldata.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000005_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -347,8 +347,8 @@ end
 
       @timed_testset "reinterpolate with nonuniform data" begin
         # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir)
-        outfilename = "solution_000005.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -363,8 +363,8 @@ end
 
       @timed_testset "do not reinterpolate with nonuniform data" begin
         # Create and test output without reinterpolation on LGL nodes
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir, reinterpolate=false)
-        outfilename = "solution_000005.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -380,8 +380,8 @@ end
       @timed_testset "do not reinterpolate with uniform data" begin
         # Create and test output without reinterpolation on uniform nodes
         # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000005.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-        outfilename = "solution_000005.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000005.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000005.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -481,6 +481,48 @@ end
       end
     end
   end
+
+  if VERSION >= v"1.8"
+    # Julia v1.7 heavily downgrades Trixi.jl. Subcell limiting is not yet supported.
+    # Therefore, only perform tests with Julia v1.8 or newer.
+    @testset "Subcell limiting coefficients" begin
+      isdir(outdir) && rm(outdir, recursive=true)
+      run_trixi(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
+                maxiters=10, initial_refinement_level=4)
+
+      @timed_testset "without reinterpolation" begin
+        # Create and test output without reinterpolation
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000010.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-treemesh-shockcapturing-subcell-no-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/treemesh/dgsem_sedov_subcell_no_interp_10.vtu"
+        ref_file = get_test_reference_file("dgsem_sedov_subcell_no_interp_10.vtu", remote_filename)
+        compare_point_data(out_file, ref_file)
+      end
+
+      @timed_testset "with reinterpolation" begin
+        # Create and test output without reinterpolation
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000010.h5"), output_directory=outdir, reinterpolate=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000010.vtu"
+        out_file = joinpath(outdir, outfilename)
+
+        # save output file to `artifacts` to facilitate debugging of failing tests
+        testname = "2d-treemesh-shockcapturing-subcell-reinterp"
+        cp(out_file, joinpath(artifacts_dir, testname * "-" * outfilename), force=true)
+
+        # remote file path is actually a URL so it always has the same path structure
+        remote_filename = "2d/treemesh/dgsem_sedov_subcell_interp_10.vtu"
+        ref_file = get_test_reference_file("dgsem_sedov_subcell_interp_10.vtu", remote_filename)
+        compare_cell_data(out_file, ref_file)
+      end
+    end
+  end
 end
 
 # Clean up afterwards: delete Trixi output directory and reference file directory
diff --git a/test/test_3d.jl b/test/test_3d.jl
index 5a60e07..dfe93e6 100644
--- a/test/test_3d.jl
+++ b/test/test_3d.jl
@@ -22,8 +22,8 @@ end
 
     @timed_testset "mesh data" begin
       # create the output file to be tested
-      @test_nowarn trixi2vtk(joinpath(outdir, "mesh_000004.h5"), output_directory=outdir)
-      outfilename = "mesh_000004_celldata.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "mesh_" * LEADING_ZEROS * "000004.h5"), output_directory=outdir)
+      outfilename = "mesh_" * LEADING_ZEROS * "000004_celldata.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -38,8 +38,8 @@ end
 
     @timed_testset "solution celldata" begin
       # create the output file to be tested
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir)
-      outfilename = "solution_000004_celldata.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000004.h5"), output_directory=outdir)
+      outfilename = "solution_" * LEADING_ZEROS * "000004_celldata.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -54,8 +54,8 @@ end
 
     @timed_testset "reinterpolate with nonuniform data" begin
       # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir)
-      outfilename = "solution_000004.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000004.h5"), output_directory=outdir)
+      outfilename = "solution_" * LEADING_ZEROS * "000004.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -70,8 +70,8 @@ end
 
     @timed_testset "do not reinterpolate with nonuniform data" begin
       # Create and test output without reinterpolation on LGL nodes
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir, reinterpolate=false)
-      outfilename = "solution_000004.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000004.h5"), output_directory=outdir, reinterpolate=false)
+      outfilename = "solution_" * LEADING_ZEROS * "000004.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -87,8 +87,8 @@ end
     @timed_testset "do not reinterpolate with uniform data" begin
       # Create and test output without reinterpolation on uniform nodes
       # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-      @test_nowarn trixi2vtk(joinpath(outdir, "solution_000004.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-      outfilename = "solution_000004.vtu"
+      @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000004.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+      outfilename = "solution_" * LEADING_ZEROS * "000004.vtu"
       out_file = joinpath(outdir, outfilename)
 
       # save output file to `artifacts` to facilitate debugging of failing tests
@@ -127,8 +127,8 @@ end
 
       @timed_testset "solution celldata" begin
         # create the output file to be tested
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir)
-        outfilename = "solution_000001_celldata.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000001_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -143,8 +143,8 @@ end
 
       @timed_testset "reinterpolate with nonuniform data" begin
         # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -159,8 +159,8 @@ end
 
       @timed_testset "do not reinterpolate with nonuniform data" begin
         # Create and test output without reinterpolation on LGL nodes
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -176,8 +176,8 @@ end
       @timed_testset "do not reinterpolate with uniform data" begin
         # Create and test output without reinterpolation on uniform nodes
         # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-        outfilename = "solution_000001.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000001.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000001.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -213,8 +213,8 @@ end
 
       @timed_testset "solution celldata" begin
         # create the output file to be tested
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir)
-        outfilename = "solution_000002_celldata.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000002_celldata.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -229,8 +229,8 @@ end
 
       @timed_testset "reinterpolate with nonuniform data" begin
         # Create and test output with reinterpolation (default options: `reinterpolate=true, data_is_uniform=false`)
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir)
-        outfilename = "solution_000002.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -245,8 +245,8 @@ end
 
       @timed_testset "do not reinterpolate with nonuniform data" begin
         # Create and test output without reinterpolation on LGL nodes
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir, reinterpolate=false)
-        outfilename = "solution_000002.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir, reinterpolate=false)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
@@ -262,8 +262,8 @@ end
       @timed_testset "do not reinterpolate with uniform data" begin
         # Create and test output without reinterpolation on uniform nodes
         # OBS! This is a dummy test just to exercise code. The resulting plot will look weird.
-        @test_nowarn trixi2vtk(joinpath(outdir, "solution_000002.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
-        outfilename = "solution_000002.vtu"
+        @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000002.h5"), output_directory=outdir, reinterpolate=false, data_is_uniform=true)
+        outfilename = "solution_" * LEADING_ZEROS * "000002.vtu"
         out_file = joinpath(outdir, outfilename)
 
         # save output file to `artifacts` to facilitate debugging of failing tests
diff --git a/test/test_manual.jl b/test/test_manual.jl
index a051ee4..79f1d03 100644
--- a/test/test_manual.jl
+++ b/test/test_manual.jl
@@ -29,7 +29,7 @@ isdir(outdir) && rm(outdir, recursive=true)
     end
 
     @testset "unsupported file format" begin
-      @test_throws ErrorException trixi2vtk(joinpath(outdir, "solution_000000.h5");
+      @test_throws ErrorException trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000000.h5");
                                             output_directory=outdir,
                                             format=:does_not_exist)
     end
@@ -40,9 +40,9 @@ isdir(outdir) && rm(outdir, recursive=true)
   end
 
   @testset "trixi2vtk set number of output nodes" begin
-    @test_nowarn trixi2vtk(joinpath(outdir, "solution_000000.h5"); nvisnodes=0)
+    @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000000.h5"); nvisnodes=0)
 
-    @test_nowarn trixi2vtk(joinpath(outdir, "solution_000000.h5"); nvisnodes=5)
+    @test_nowarn trixi2vtk(joinpath(outdir, "solution_" * LEADING_ZEROS * "000000.h5"); nvisnodes=5)
   end
 
   @timed_testset "pvd_filenames" begin
diff --git a/test/test_trixi2vtk.jl b/test/test_trixi2vtk.jl
index 15e59ad..b940dfb 100644
--- a/test/test_trixi2vtk.jl
+++ b/test/test_trixi2vtk.jl
@@ -5,6 +5,21 @@ using Trixi2Vtk
 using ReadVTK
 
 
+# Get the version of Trixi.jl we are testing since the output file name
+# format changed in v0.8.0
+import Pkg
+using UUIDs: UUID
+const LEADING_ZEROS = let
+  trixi_uuid = UUID("a7f1ee26-1774-49b1-8366-f1abc58fbfcb")
+  trixi_version = Pkg.dependencies()[trixi_uuid].version
+  if trixi_version > v"0.7.999999999"
+    "000"
+  else
+    ""
+  end
+end
+
+
 function run_trixi(elixir; kwargs...)
   # evaluate examples in the scope of the module they're called from
   trixi_include(@__MODULE__, elixir; kwargs...)
@@ -41,7 +56,7 @@ end
 if VERSION < v"1.8"
   const TEST_REFERENCE_COMMIT = "c0a966b06489f9b2ee3aefeb0a5c0dae733df36f"
 else
-  const TEST_REFERENCE_COMMIT = "86a43fe8dc254130345754fb512268204cf2233c"
+  const TEST_REFERENCE_COMMIT = "e51f3613ac1adfcfd2bf1d74a2756034dab0579c"
 end
 
 # Local folder to store downloaded reference files. If you change this, also adapt `../.gitignore`!