diff --git a/test/ParallelIOTest.cpp b/test/ParallelIOTest.cpp
index b6621ea8f1..cce13ba59d 100644
--- a/test/ParallelIOTest.cpp
+++ b/test/ParallelIOTest.cpp
@@ -1766,4 +1766,160 @@ TEST_CASE("unavailable_backend", "[core][parallel]")
     }
 #endif
 }
+
+void joined_dim(std::string const &ext)
+{
+    using type = float;
+    using patchType = uint64_t;
+    constexpr size_t patches_per_rank = 5;
+    constexpr size_t length_of_patch = 10;
+
+    int size{-1};
+    int rank{-1};
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    {
+        Series s(
+            "../samples/joinedDimParallel." + ext,
+            Access::CREATE,
+            MPI_COMM_WORLD);
+        std::vector<UniquePtrWithLambda<type>> writeFrom(patches_per_rank);
+
+        auto it = s.writeIterations()[100];
+
+        Dataset numParticlesDS(
+            determineDatatype<patchType>(), {Dataset::JOINED_DIMENSION});
+        auto numParticles =
+            it.particles["e"]
+                .particlePatches["numParticles"][RecordComponent::SCALAR];
+        auto numParticlesOffset =
+            it.particles["e"]
+                .particlePatches["numParticlesOffset"][RecordComponent::SCALAR];
+        numParticles.resetDataset(numParticlesDS);
+        numParticlesOffset.resetDataset(numParticlesDS);
+
+        auto patchOffset = it.particles["e"].particlePatches["offset"]["x"];
+        auto patchExtent = it.particles["e"].particlePatches["extent"]["x"];
+        Dataset particlePatchesDS(
+            determineDatatype<float>(), {Dataset::JOINED_DIMENSION});
+        patchOffset.resetDataset(particlePatchesDS);
+        patchExtent.resetDataset(particlePatchesDS);
+
+        float start_value = rank * patches_per_rank * length_of_patch;
+        for (size_t i = 0; i < 5; ++i)
+        {
+            writeFrom[i] = UniquePtrWithLambda<type>(
+                new type[length_of_patch],
+                [](auto const *ptr) { delete[] ptr; });
+            std::iota(
+                writeFrom[i].get(),
+                writeFrom[i].get() + 10,
+                start_value + length_of_patch * i);
+            patchOffset.store<type>(start_value + length_of_patch * i);
+        }
+
+        auto epx = it.particles["e"]["position"]["x"];
+        Dataset ds(determineDatatype<type>(), {Dataset::JOINED_DIMENSION});
+        epx.resetDataset(ds);
+
+        size_t counter = 0;
+        for (auto &chunk : writeFrom)
+        {
+            epx.storeChunk(std::move(chunk), {}, {length_of_patch});
+            numParticles.store<patchType>(length_of_patch);
+            /*
+             * For the sake of the test case, we know that the
+             * numParticlesOffset has this value. In general, the purpose of the
+             * joined array is that we don't need to know these values, so the
+             * specification of particle patches is somewhat difficult.
+             */
+            numParticlesOffset.store<patchType>(
+                start_value + counter++ * length_of_patch);
+            patchExtent.store<type>(10);
+        }
+        writeFrom.clear();
+        it.close();
+        s.close();
+    }
+
+    {
+        Series s(
+            "../samples/joinedDimParallel." + ext,
+            Access::READ_ONLY,
+            MPI_COMM_WORLD);
+        auto it = s.iterations[100];
+        auto e = it.particles["e"];
+
+        auto particleData = e["position"]["x"].loadChunk<type>();
+        auto numParticles =
+            e.particlePatches["numParticles"][RecordComponent::SCALAR]
+                .load<patchType>();
+        auto numParticlesOffset =
+            e.particlePatches["numParticlesOffset"][RecordComponent::SCALAR]
+                .load<patchType>();
+        auto patchOffset = e.particlePatches["offset"]["x"].load<type>();
+        auto patchExtent = e.particlePatches["extent"]["x"].load<type>();
+
+        it.close();
+
+        // check validity of particle patches
+        auto numPatches =
+            e.particlePatches["numParticlesOffset"][RecordComponent::SCALAR]
+                .getExtent()[0];
+        REQUIRE(
+            e.particlePatches["numParticles"][RecordComponent::SCALAR]
+                .getExtent()[0] == numPatches);
+        for (size_t i = 0; i < numPatches; ++i)
+        {
+            for (size_t j = 0; j < numParticles.get()[i]; ++j)
+            {
+                REQUIRE(
+                    patchOffset.get()[i] <=
+                    particleData.get()[numParticlesOffset.get()[i] + j]);
+                REQUIRE(
+                    particleData.get()[numParticlesOffset.get()[i] + j] <
+                    patchOffset.get()[i] + patchExtent.get()[i]);
+            }
+        }
+
+        /*
+         * Check that joined array joins early writes before later writes from
+         * the same rank
+         */
+        for (size_t i = 0; i < size * length_of_patch * patches_per_rank; ++i)
+        {
+            REQUIRE(float(i) == particleData.get()[i]);
+        }
+        for (size_t i = 0; i < size * patches_per_rank; ++i)
+        {
+            REQUIRE(length_of_patch * i == numParticlesOffset.get()[i]);
+            REQUIRE(type(length_of_patch * i) == patchOffset.get()[i]);
+        }
+    }
+}
+
+TEST_CASE("joined_dim", "[parallel]")
+{
+#if 100000000 * ADIOS2_VERSION_MAJOR + 1000000 * ADIOS2_VERSION_MINOR +        \
+        10000 * ADIOS2_VERSION_PATCH + 100 * ADIOS2_VERSION_TWEAK >=           \
+    209000000
+    constexpr char const *supportsJoinedDims[] = {"bp", "bp4", "bp5"};
+#else
+    // no zero-size arrays
+    std::vector<char const *> supportsJoinedDims;
+#endif
+    for (auto const &t : testedFileExtensions())
+    {
+        for (auto const supported : supportsJoinedDims)
+        {
+            if (t == supported)
+            {
+                joined_dim(t);
+                break;
+            }
+        }
+    }
+}
+
 #endif // openPMD_HAVE_ADIOS2 && openPMD_HAVE_MPI
diff --git a/test/SerialIOTest.cpp b/test/SerialIOTest.cpp
index ce37e4f352..f0b6f9557e 100644
--- a/test/SerialIOTest.cpp
+++ b/test/SerialIOTest.cpp
@@ -183,13 +183,13 @@ TEST_CASE("adios2_char_portability", "[serial][adios2]")
         REQUIRE(unspecifiedVectorAttribute.dtype == Datatype::VEC_STRING);
         std::vector<std::string> desiredVector{"abc", "def", "ghi"};
         REQUIRE(
-            signedVectorAttribute.get<std::vector<std::string> >() ==
+            signedVectorAttribute.get<std::vector<std::string>>() ==
             desiredVector);
         REQUIRE(
-            unsignedVectorAttribute.get<std::vector<std::string> >() ==
+            unsignedVectorAttribute.get<std::vector<std::string>>() ==
             desiredVector);
         REQUIRE(
-            unspecifiedVectorAttribute.get<std::vector<std::string> >() ==
+            unspecifiedVectorAttribute.get<std::vector<std::string>>() ==
             desiredVector);
 
         auto signedCharAttribute = read.getAttribute("signedChar");
@@ -998,7 +998,7 @@ inline void constant_scalar(std::string file_ending)
             s.iterations[1]
                 .meshes["rho"][MeshRecordComponent::SCALAR]
                 .getAttribute("shape")
-                .get<std::vector<uint64_t> >() == Extent{1, 2, 3});
+                .get<std::vector<uint64_t>>() == Extent{1, 2, 3});
         REQUIRE(s.iterations[1]
                     .meshes["rho"][MeshRecordComponent::SCALAR]
                     .containsAttribute("value"));
@@ -1020,7 +1020,7 @@ inline void constant_scalar(std::string file_ending)
             s.iterations[1]
                 .meshes["E"]["x"]
                 .getAttribute("shape")
-                .get<std::vector<uint64_t> >() == Extent{1, 2, 3});
+                .get<std::vector<uint64_t>>() == Extent{1, 2, 3});
         REQUIRE(s.iterations[1].meshes["E"]["x"].containsAttribute("value"));
         REQUIRE(
             s.iterations[1]
@@ -1048,7 +1048,7 @@ inline void constant_scalar(std::string file_ending)
             s.iterations[1]
                 .particles["e"]["position"][RecordComponent::SCALAR]
                 .getAttribute("shape")
-                .get<std::vector<uint64_t> >() == Extent{3, 2, 1});
+                .get<std::vector<uint64_t>>() == Extent{3, 2, 1});
         REQUIRE(s.iterations[1]
                     .particles["e"]["position"][RecordComponent::SCALAR]
                     .containsAttribute("value"));
@@ -1072,7 +1072,7 @@ inline void constant_scalar(std::string file_ending)
             s.iterations[1]
                 .particles["e"]["positionOffset"][RecordComponent::SCALAR]
                 .getAttribute("shape")
-                .get<std::vector<uint64_t> >() == Extent{3, 2, 1});
+                .get<std::vector<uint64_t>>() == Extent{3, 2, 1});
         REQUIRE(s.iterations[1]
                     .particles["e"]["positionOffset"][RecordComponent::SCALAR]
                     .containsAttribute("value"));
@@ -1094,7 +1094,7 @@ inline void constant_scalar(std::string file_ending)
             s.iterations[1]
                 .particles["e"]["velocity"]["x"]
                 .getAttribute("shape")
-                .get<std::vector<uint64_t> >() == Extent{3, 2, 1});
+                .get<std::vector<uint64_t>>() == Extent{3, 2, 1});
         REQUIRE(
             s.iterations[1].particles["e"]["velocity"]["x"].containsAttribute(
                 "value"));
@@ -1446,55 +1446,55 @@ inline void dtype_test(const std::string &backend)
         REQUIRE(s.getAttribute("emptyString").get<std::string>().empty());
     }
     REQUIRE(
-        s.getAttribute("vecChar").get<std::vector<char> >() ==
+        s.getAttribute("vecChar").get<std::vector<char>>() ==
         std::vector<char>({'c', 'h', 'a', 'r'}));
     REQUIRE(
-        s.getAttribute("vecInt16").get<std::vector<int16_t> >() ==
+        s.getAttribute("vecInt16").get<std::vector<int16_t>>() ==
         std::vector<int16_t>({32766, 32767}));
     REQUIRE(
-        s.getAttribute("vecInt32").get<std::vector<int32_t> >() ==
+        s.getAttribute("vecInt32").get<std::vector<int32_t>>() ==
         std::vector<int32_t>({2147483646, 2147483647}));
     REQUIRE(
-        s.getAttribute("vecInt64").get<std::vector<int64_t> >() ==
+        s.getAttribute("vecInt64").get<std::vector<int64_t>>() ==
         std::vector<int64_t>({9223372036854775806, 9223372036854775807}));
     REQUIRE(
-        s.getAttribute("vecUchar").get<std::vector<unsigned char> >() ==
+        s.getAttribute("vecUchar").get<std::vector<unsigned char>>() ==
         std::vector<unsigned char>({'u', 'c', 'h', 'a', 'r'}));
     REQUIRE(
-        s.getAttribute("vecSchar").get<std::vector<signed char> >() ==
+        s.getAttribute("vecSchar").get<std::vector<signed char>>() ==
         std::vector<signed char>({'s', 'c', 'h', 'a', 'r'}));
     REQUIRE(
-        s.getAttribute("vecUint16").get<std::vector<uint16_t> >() ==
+        s.getAttribute("vecUint16").get<std::vector<uint16_t>>() ==
         std::vector<uint16_t>({65534u, 65535u}));
     REQUIRE(
-        s.getAttribute("vecUint32").get<std::vector<uint32_t> >() ==
+        s.getAttribute("vecUint32").get<std::vector<uint32_t>>() ==
         std::vector<uint32_t>({4294967294u, 4294967295u}));
     REQUIRE(
-        s.getAttribute("vecUint64").get<std::vector<uint64_t> >() ==
+        s.getAttribute("vecUint64").get<std::vector<uint64_t>>() ==
         std::vector<uint64_t>({18446744073709551614u, 18446744073709551615u}));
     REQUIRE(
-        s.getAttribute("vecFloat").get<std::vector<float> >() ==
+        s.getAttribute("vecFloat").get<std::vector<float>>() ==
         std::vector<float>({0.f, 3.40282e+38f}));
     REQUIRE(
-        s.getAttribute("vecDouble").get<std::vector<double> >() ==
+        s.getAttribute("vecDouble").get<std::vector<double>>() ==
         std::vector<double>({0., 1.79769e+308}));
     if (test_long_double)
     {
         REQUIRE(
-            s.getAttribute("vecLongdouble").get<std::vector<long double> >() ==
+            s.getAttribute("vecLongdouble").get<std::vector<long double>>() ==
             std::vector<long double>(
                 {0.L, std::numeric_limits<long double>::max()}));
     }
     REQUIRE(
-        s.getAttribute("vecString").get<std::vector<std::string> >() ==
+        s.getAttribute("vecString").get<std::vector<std::string>>() ==
         std::vector<std::string>({"vector", "of", "strings"}));
     if (!adios1)
     {
         REQUIRE(
-            s.getAttribute("vecEmptyString").get<std::vector<std::string> >() ==
+            s.getAttribute("vecEmptyString").get<std::vector<std::string>>() ==
             std::vector<std::string>({"", "", ""}));
         REQUIRE(
-            s.getAttribute("vecMixedString").get<std::vector<std::string> >() ==
+            s.getAttribute("vecMixedString").get<std::vector<std::string>>() ==
             std::vector<std::string>({"hi", "", "ho"}));
     }
     REQUIRE(s.getAttribute("bool").get<bool>() == true);
@@ -1691,7 +1691,7 @@ void test_complex(const std::string &backend)
                 "longDoublesYouSay", std::complex<long double>(5.5, -4.55));
 
         auto Cflt = o.iterations[0].meshes["Cflt"][RecordComponent::SCALAR];
-        std::vector<std::complex<float> > cfloats(3);
+        std::vector<std::complex<float>> cfloats(3);
         cfloats.at(0) = {1., 2.};
         cfloats.at(1) = {-3., 4.};
         cfloats.at(2) = {5., -6.};
@@ -1699,14 +1699,14 @@ void test_complex(const std::string &backend)
         Cflt.storeChunk(cfloats, {0});
 
         auto Cdbl = o.iterations[0].meshes["Cdbl"][RecordComponent::SCALAR];
-        std::vector<std::complex<double> > cdoubles(3);
+        std::vector<std::complex<double>> cdoubles(3);
         cdoubles.at(0) = {2., 1.};
         cdoubles.at(1) = {-4., 3.};
         cdoubles.at(2) = {6., -5.};
         Cdbl.resetDataset(Dataset(Datatype::CDOUBLE, {cdoubles.size()}));
         Cdbl.storeChunk(cdoubles, {0});
 
-        std::vector<std::complex<long double> > cldoubles(3);
+        std::vector<std::complex<long double>> cldoubles(3);
         if (o.backend() != "ADIOS2" && o.backend() != "ADIOS1" &&
             o.backend() != "MPI_ADIOS1")
         {
@@ -1727,26 +1727,26 @@ void test_complex(const std::string &backend)
         Series i = Series(
             "../samples/serial_write_complex." + backend, Access::READ_ONLY);
         REQUIRE(
-            i.getAttribute("lifeIsComplex").get<std::complex<double> >() ==
+            i.getAttribute("lifeIsComplex").get<std::complex<double>>() ==
             std::complex<double>(4.56, 7.89));
         REQUIRE(
-            i.getAttribute("butComplexFloats").get<std::complex<float> >() ==
+            i.getAttribute("butComplexFloats").get<std::complex<float>>() ==
             std::complex<float>(42.3, -99.3));
         if (i.backend() != "ADIOS2" && i.backend() != "ADIOS1" &&
             i.backend() != "MPI_ADIOS1")
         {
             REQUIRE(
                 i.getAttribute("longDoublesYouSay")
-                    .get<std::complex<long double> >() ==
+                    .get<std::complex<long double>>() ==
                 std::complex<long double>(5.5, -4.55));
         }
 
         auto rcflt = i.iterations[0]
                          .meshes["Cflt"][RecordComponent::SCALAR]
-                         .loadChunk<std::complex<float> >();
+                         .loadChunk<std::complex<float>>();
         auto rcdbl = i.iterations[0]
                          .meshes["Cdbl"][RecordComponent::SCALAR]
-                         .loadChunk<std::complex<double> >();
+                         .loadChunk<std::complex<double>>();
         i.flush();
 
         REQUIRE(rcflt.get()[1] == std::complex<float>(-3., 4.));
@@ -1757,7 +1757,7 @@ void test_complex(const std::string &backend)
         {
             auto rcldbl = i.iterations[0]
                               .meshes["Cldbl"][RecordComponent::SCALAR]
-                              .loadChunk<std::complex<long double> >();
+                              .loadChunk<std::complex<long double>>();
             i.flush();
             REQUIRE(rcldbl.get()[2] == std::complex<long double>(7., -6.));
         }
@@ -4907,7 +4907,7 @@ void bp4_steps(
         auto E_x = E["x"];
         REQUIRE(
             E.getAttribute("vector_of_string")
-                .get<std::vector<std::string> >() ==
+                .get<std::vector<std::string>>() ==
             std::vector<std::string>{"vector", "of", "string"});
         REQUIRE(E_x.getDimensionality() == 1);
         REQUIRE(E_x.getExtent()[0] == 10);
@@ -5133,7 +5133,7 @@ struct AreEqual
 };
 
 template <typename T>
-struct AreEqual<std::vector<T> >
+struct AreEqual<std::vector<T>>
 {
     static bool areEqual(std::vector<T> v1, std::vector<T> v2)
     {
@@ -7215,3 +7215,147 @@ TEST_CASE("groupbased_read_write", "[serial]")
         }
     }
 }
+
+void joined_dim(std::string const &ext)
+{
+    using type = float;
+    using patchType = uint64_t;
+    constexpr size_t patches_per_rank = 5;
+    constexpr size_t length_of_patch = 10;
+
+    {
+        Series s("../samples/joinedDimParallel." + ext, Access::CREATE);
+        std::vector<UniquePtrWithLambda<type>> writeFrom(patches_per_rank);
+
+        auto it = s.writeIterations()[100];
+
+        Dataset numParticlesDS(
+            determineDatatype<patchType>(), {Dataset::JOINED_DIMENSION});
+        auto numParticles =
+            it.particles["e"]
+                .particlePatches["numParticles"][RecordComponent::SCALAR];
+        auto numParticlesOffset =
+            it.particles["e"]
+                .particlePatches["numParticlesOffset"][RecordComponent::SCALAR];
+        numParticles.resetDataset(numParticlesDS);
+        numParticlesOffset.resetDataset(numParticlesDS);
+
+        auto patchOffset = it.particles["e"].particlePatches["offset"]["x"];
+        auto patchExtent = it.particles["e"].particlePatches["extent"]["x"];
+        Dataset particlePatchesDS(
+            determineDatatype<float>(), {Dataset::JOINED_DIMENSION});
+        patchOffset.resetDataset(particlePatchesDS);
+        patchExtent.resetDataset(particlePatchesDS);
+
+        for (size_t i = 0; i < 5; ++i)
+        {
+            writeFrom[i] = UniquePtrWithLambda<type>(
+                new type[length_of_patch],
+                [](auto const *ptr) { delete[] ptr; });
+            std::iota(
+                writeFrom[i].get(),
+                writeFrom[i].get() + 10,
+                length_of_patch * i);
+            patchOffset.store<type>(length_of_patch * i);
+        }
+
+        auto epx = it.particles["e"]["position"]["x"];
+        Dataset ds(determineDatatype<type>(), {Dataset::JOINED_DIMENSION});
+        epx.resetDataset(ds);
+
+        size_t counter = 0;
+        for (auto &chunk : writeFrom)
+        {
+            epx.storeChunk(std::move(chunk), {}, {length_of_patch});
+            numParticles.store<patchType>(length_of_patch);
+            /*
+             * For the sake of the test case, we know that the
+             * numParticlesOffset has this value. In general, the purpose of the
+             * joined array is that we don't need to know these values, so the
+             * specification of particle patches is somewhat difficult.
+             */
+            numParticlesOffset.store<patchType>(counter++ * length_of_patch);
+            patchExtent.store<type>(10);
+        }
+        writeFrom.clear();
+        it.close();
+        s.close();
+    }
+
+    {
+        Series s("../samples/joinedDimParallel." + ext, Access::READ_ONLY);
+        auto it = s.iterations[100];
+        auto e = it.particles["e"];
+
+        auto particleData = e["position"]["x"].loadChunk<type>();
+        auto numParticles =
+            e.particlePatches["numParticles"][RecordComponent::SCALAR]
+                .load<patchType>();
+        auto numParticlesOffset =
+            e.particlePatches["numParticlesOffset"][RecordComponent::SCALAR]
+                .load<patchType>();
+        auto patchOffset = e.particlePatches["offset"]["x"].load<type>();
+        auto patchExtent = e.particlePatches["extent"]["x"].load<type>();
+
+        it.close();
+
+        // check validity of particle patches
+        auto numPatches =
+            e.particlePatches["numParticlesOffset"][RecordComponent::SCALAR]
+                .getExtent()[0];
+        REQUIRE(
+            e.particlePatches["numParticles"][RecordComponent::SCALAR]
+                .getExtent()[0] == numPatches);
+        for (size_t i = 0; i < numPatches; ++i)
+        {
+            for (size_t j = 0; j < numParticles.get()[i]; ++j)
+            {
+                REQUIRE(
+                    patchOffset.get()[i] <=
+                    particleData.get()[numParticlesOffset.get()[i] + j]);
+                REQUIRE(
+                    particleData.get()[numParticlesOffset.get()[i] + j] <
+                    patchOffset.get()[i] + patchExtent.get()[i]);
+            }
+        }
+
+        /*
+         * Check that:
+         * 1. Joined array joins writes from lower ranks before higher ranks
+         * 2. Joined array joins early writes before later writes from the same
+         *    rank
+         */
+        for (size_t i = 0; i < length_of_patch * patches_per_rank; ++i)
+        {
+            REQUIRE(float(i) == particleData.get()[i]);
+        }
+        for (size_t i = 0; i < patches_per_rank; ++i)
+        {
+            REQUIRE(length_of_patch * i == numParticlesOffset.get()[i]);
+            REQUIRE(type(length_of_patch * i) == patchOffset.get()[i]);
+        }
+    }
+}
+
+TEST_CASE("joined_dim", "[serial]")
+{
+#if 100000000 * ADIOS2_VERSION_MAJOR + 1000000 * ADIOS2_VERSION_MINOR +        \
+        10000 * ADIOS2_VERSION_PATCH + 100 * ADIOS2_VERSION_TWEAK >=           \
+    209000000
+    constexpr char const *supportsJoinedDims[] = {"bp", "bp4", "bp5"};
+#else
+    // no zero-size arrays
+    std::vector<char const *> supportsJoinedDims;
+#endif
+    for (auto const &t : testedFileExtensions())
+    {
+        for (auto const supported : supportsJoinedDims)
+        {
+            if (t == supported)
+            {
+                joined_dim(t);
+                break;
+            }
+        }
+    }
+}