diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1a1b21d87c..64485c3b07 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -243,7 +243,7 @@ endif()
 
 set(ADIOS2_CONFIG_OPTS
     DataMan DataSpaces HDF5 HDF5_VOL MHS SST Fortran MPI Python Blosc2 BZip2
-    LIBPRESSIO MGARD PNG SZ ZFP DAOS IME O_DIRECT Sodium Catalyst SysVShMem UCX
+    LIBPRESSIO MGARD MGARD_MDR PNG SZ ZFP DAOS IME O_DIRECT Sodium Catalyst SysVShMem UCX
     ZeroMQ Profiling Endian_Reverse Derived_Variable AWSSDK GPU_Support CUDA Kokkos
     Kokkos_CUDA Kokkos_HIP Kokkos_SYCL
 )
diff --git a/bindings/CXX11/adios2/cxx11/Variable.cpp b/bindings/CXX11/adios2/cxx11/Variable.cpp
index cc1944fae9..d7d47c3428 100644
--- a/bindings/CXX11/adios2/cxx11/Variable.cpp
+++ b/bindings/CXX11/adios2/cxx11/Variable.cpp
@@ -78,6 +78,13 @@ namespace adios2
     }                                                                                              \
                                                                                                    \
     template <>                                                                                    \
+    void Variable<T>::SetAccuracy(const adios2::Accuracy &a)                                       \
+    {                                                                                              \
+        helper::CheckForNullptr(m_Variable, "in call to Variable<T>::SetAccuracy");                \
+        m_Variable->SetAccuracy(a);                                                                \
+    }                                                                                              \
+                                                                                                   \
+    template <>                                                                                    \
     size_t Variable<T>::SelectionSize() const                                                      \
     {                                                                                              \
         helper::CheckForNullptr(m_Variable, "in call to Variable<T>::SelectionSize");              \
@@ -219,6 +226,13 @@ namespace adios2
     }                                                                                              \
                                                                                                    \
     template <>                                                                                    \
+    adios2::Accuracy Variable<T>::GetAccuracy()                                                    \
+    {                                                                                              \
+        helper::CheckForNullptr(m_Variable, "in call to Variable<T>::GetAccuracy");                \
+        return m_Variable->GetAccuracy();                                                          \
+    }                                                                                              \
+                                                                                                   \
+    template <>                                                                                    \
     std::vector<std::vector<typename Variable<T>::Info>> Variable<T>::AllStepsBlocksInfo()         \
     {                                                                                              \
         return DoAllStepsBlocksInfo();                                                             \
diff --git a/bindings/CXX11/adios2/cxx11/Variable.h b/bindings/CXX11/adios2/cxx11/Variable.h
index 398d001667..4f8ea8f071 100644
--- a/bindings/CXX11/adios2/cxx11/Variable.h
+++ b/bindings/CXX11/adios2/cxx11/Variable.h
@@ -209,6 +209,12 @@ class Variable
      */
     void SetStepSelection(const adios2::Box<size_t> &stepSelection);
 
+    /**
+     * Sets the requested accuracy for the next read operation.
+     * The actual accuracy after the read is provided in GetAccuracy()
+     */
+    void SetAccuracy(const adios2::Accuracy &a);
+
     /**
      * Returns the number of elements required for pre-allocation based on
      * current count and stepsCount
@@ -335,6 +341,13 @@ class Variable
      */
     T Max(const size_t step = adios2::DefaultSizeT) const;
 
+    /**
+     * Get the provided accuracy for the last read operation.
+     * Most operations provide data as it was written, meaning that
+     * error is reported as 0.0
+     */
+    adios2::Accuracy GetAccuracy();
+
     /** Contains block information for a particular Variable<T> */
     struct Info
     {
diff --git a/cmake/DetectOptions.cmake b/cmake/DetectOptions.cmake
index d66e096bda..4dd7f5522a 100644
--- a/cmake/DetectOptions.cmake
+++ b/cmake/DetectOptions.cmake
@@ -175,6 +175,7 @@ if(TARGET LibPressio::libpressio)
 endif()
 
 # MGARD
+set(ADIOS2_HAVE_MGARD_MDR FALSE)
 if(ADIOS2_USE_MGARD STREQUAL AUTO)
   find_package(mgard CONFIG)
 elseif(ADIOS2_USE_MGARD)
@@ -182,6 +183,9 @@ elseif(ADIOS2_USE_MGARD)
 endif()
 if(mgard_FOUND)
   set(ADIOS2_HAVE_MGARD TRUE)
+  if(MGARD_ENABLE_MDR)
+     set(ADIOS2_HAVE_MGARD_MDR TRUE)
+ endif()
 endif()
 
 # PNG
diff --git a/examples/basics/globalArrayND/globalArrayNDWrite.cpp b/examples/basics/globalArrayND/globalArrayNDWrite.cpp
index b1509b7b54..7658a58843 100644
--- a/examples/basics/globalArrayND/globalArrayNDWrite.cpp
+++ b/examples/basics/globalArrayND/globalArrayNDWrite.cpp
@@ -49,7 +49,7 @@ int main(int argc, char *argv[])
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
     MPI_Comm_size(MPI_COMM_WORLD, &nproc);
 #endif
-    const int NSTEPS = 5;
+    const int NSTEPS = 1;
 
 #if ADIOS2_USE_MPI
     adios2::ADIOS adios(MPI_COMM_WORLD);
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
 #endif
 
     // Application variables for output
-    const unsigned int Nx = 10;
+    const unsigned int Nx = 100000;
     // Global 2D array, size of nproc x Nx, with 1D decomposition
     // Each process writes one "row" of the 2D matrix.
     std::vector<double> row(Nx);
@@ -84,6 +84,9 @@ int main(int argc, char *argv[])
 
         io.DefineAttribute<int>("nsteps", NSTEPS);
 
+        adios2::Operator op = adios.DefineOperator("mdr", "mdr");
+        varGlobalArray.AddOperation(op, {{"accuracy", std::to_string(0.1)}});
+
         // Open file. "w" means we overwrite any existing file on disk,
         // but Advance() will append steps to the same file.
         adios2::Engine writer = io.Open("globalArray.bp", adios2::Mode::Write);
diff --git a/source/adios2/CMakeLists.txt b/source/adios2/CMakeLists.txt
index 355db2992d..637d8ea361 100644
--- a/source/adios2/CMakeLists.txt
+++ b/source/adios2/CMakeLists.txt
@@ -323,6 +323,9 @@ endif()
 if(ADIOS2_HAVE_MGARD)
   target_sources(adios2_core PRIVATE operator/compress/CompressMGARDPlus.cpp)
   target_sources(adios2_core PRIVATE operator/compress/CompressMGARD.cpp)
+  if(ADIOS2_HAVE_MGARD_MDR)
+      target_sources(adios2_core PRIVATE operator/refactor/RefactorMDR.cpp)
+  endif()
   target_link_libraries(adios2_core PRIVATE mgard::mgard)
 endif()
 
diff --git a/source/adios2/common/ADIOSTypes.h b/source/adios2/common/ADIOSTypes.h
index 01423e781f..bfc0ef6df0 100644
--- a/source/adios2/common/ADIOSTypes.h
+++ b/source/adios2/common/ADIOSTypes.h
@@ -310,6 +310,22 @@ using Box = std::pair<T, T>;
 template <typename T, typename Enable = void>
 struct TypeInfo;
 
+/** Data accuracy **/
+
+/* Error. Accuracy can be requested for reading data.
+   norm: 0.0 = L2, inf() = Linf
+   relative: true = relative error, false = absolute error
+ */
+struct Accuracy
+{
+    double error;
+    double norm;
+    bool relative;
+};
+
+constexpr double L2_norm = 0.0;
+constexpr double Linf_norm = std::numeric_limits<double>::infinity();
+
 /**
  *  Return the actual size in bytes of elements of the given type.  Returns -1
  * for strings.
@@ -405,6 +421,19 @@ constexpr char s[] = "s";
 }
 #endif
 
+#ifdef ADIOS2_HAVE_MGARD_MDR
+
+constexpr char MDR[] = "mdr";
+
+namespace mdr
+{
+constexpr double DOUBLE_ROUNDING_ERROR_LIMIT = 5.0e-16;
+constexpr double FLOAT_ROUNDING_ERROR_LIMIT = 3.0e-7;
+namespace key
+{}
+}
+#endif
+
 #ifdef ADIOS2_HAVE_LIBPRESSIO
 constexpr char LossyLIBPRESSIO[] = "libpressio";
 namespace libpressio
diff --git a/source/adios2/core/Info.cpp b/source/adios2/core/Info.cpp
index 39686b8e0f..ae5dca63b3 100644
--- a/source/adios2/core/Info.cpp
+++ b/source/adios2/core/Info.cpp
@@ -72,6 +72,9 @@ static const char *const operators[] = {
     "MGARD",
     "MGARDPlus",
 #endif
+#ifdef ADIOS2_HAVE_MGARD_MDR
+    "MDR",
+#endif
 #ifdef ADIOS2_HAVE_SZ
     "SZ",
 #endif
diff --git a/source/adios2/core/Operator.cpp b/source/adios2/core/Operator.cpp
index 946c6fbada..481749e57a 100644
--- a/source/adios2/core/Operator.cpp
+++ b/source/adios2/core/Operator.cpp
@@ -11,6 +11,8 @@
 #include "Operator.h"
 #include "adios2/helper/adiosFunctions.h"
 
+#include <iostream>
+
 namespace adios2
 {
 namespace core
@@ -30,6 +32,9 @@ void Operator::SetParameter(const std::string key, const std::string value) noex
 
 Params &Operator::GetParameters() noexcept { return m_Parameters; }
 
+void Operator::SetAccuracy(const adios2::Accuracy &a) noexcept { m_AccuracyRequested = a; }
+adios2::Accuracy Operator::GetAccuracy() const noexcept { return m_AccuracyProvided; }
+
 #define declare_type(T)                                                                            \
                                                                                                    \
     void Operator::RunCallback1(const T *arg0, const std::string &arg1, const std::string &arg2,   \
@@ -96,6 +101,12 @@ Dims Operator::ConvertDims(const Dims &dimensions, const DataType type, const si
 
 size_t Operator::GetHeaderSize() const { return 0; }
 
+size_t Operator::GetEstimatedSize(const size_t ElemCount, const size_t ElemSize, const size_t ndims,
+                                  const size_t *dims) const
+{
+    return ElemCount * ElemSize + 128;
+};
+
 // PRIVATE
 void Operator::CheckCallbackType(const std::string type) const
 {
diff --git a/source/adios2/core/Operator.h b/source/adios2/core/Operator.h
index a13b596a56..d87e6788b4 100644
--- a/source/adios2/core/Operator.h
+++ b/source/adios2/core/Operator.h
@@ -37,6 +37,7 @@ class Operator
         COMPRESS_SZ = 6,
         COMPRESS_ZFP = 7,
         COMPRESS_MGARDPLUS = 8,
+        REFACTOR_MDR = 41,
         CALLBACK_SIGNATURE1 = 51,
         CALLBACK_SIGNATURE2 = 52,
         PLUGIN_INTERFACE = 53,
@@ -56,6 +57,9 @@ class Operator
 
     Params &GetParameters() noexcept;
 
+    void SetAccuracy(const adios2::Accuracy &a) noexcept;
+    adios2::Accuracy GetAccuracy() const noexcept;
+
 #define declare_type(T)                                                                            \
     virtual void RunCallback1(const T *, const std::string &, const std::string &,                 \
                               const std::string &, const size_t, const Dims &, const Dims &,       \
@@ -68,6 +72,10 @@ class Operator
 
     virtual size_t GetHeaderSize() const;
 
+    /** Give an upper bound estimate how big the transformed data could be */
+    virtual size_t GetEstimatedSize(const size_t ElemCount, const size_t ElemSize,
+                                    const size_t ndims, const size_t *dims) const;
+
     /**
      * @param dataIn
      * @param blockStart
@@ -94,6 +102,11 @@ class Operator
     /** Parameters associated with a particular Operator */
     Params m_Parameters;
 
+    /** user requested accuracy */
+    Accuracy m_AccuracyRequested = {0.0, 0.0, false};
+    /** provided accuracy */
+    Accuracy m_AccuracyProvided = {0.0, 0.0, false};
+
     /**
      * Used by lossy compressors with a limitation on complex data types or
      * dimentions Returns a adios2::Dims object that meets the requirement of a
diff --git a/source/adios2/core/VariableBase.cpp b/source/adios2/core/VariableBase.cpp
index a6430f67ba..3f76bb1f48 100644
--- a/source/adios2/core/VariableBase.cpp
+++ b/source/adios2/core/VariableBase.cpp
@@ -204,6 +204,14 @@ void VariableBase::SetMemorySelection(const Box<Dims> &memorySelection)
     m_MemoryCount = memorySelection.second;
 }
 
+void VariableBase::SetAccuracy(const adios2::Accuracy &a) noexcept
+{
+    m_AccuracyRequested = a;
+    m_AccuracyProvided = {0.0, a.norm, a.relative};
+}
+adios2::Accuracy VariableBase::GetAccuracy() const noexcept { return m_AccuracyProvided; }
+adios2::Accuracy VariableBase::GetAccuracyRequested() const noexcept { return m_AccuracyRequested; }
+
 size_t VariableBase::GetAvailableStepsStart() const { return m_AvailableStepsStart; }
 
 size_t VariableBase::GetAvailableStepsCount() const { return m_AvailableStepsCount; }
diff --git a/source/adios2/core/VariableBase.h b/source/adios2/core/VariableBase.h
index bd888a3828..b777400f40 100644
--- a/source/adios2/core/VariableBase.h
+++ b/source/adios2/core/VariableBase.h
@@ -98,6 +98,11 @@ class VariableBase
 
     Engine *m_Engine = nullptr;
 
+    /** user requested accuracy */
+    Accuracy m_AccuracyRequested = {0.0, 0.0, false};
+    /** provided accuracy */
+    Accuracy m_AccuracyProvided = {0.0, 0.0, false};
+
     /** Index to Step and blocks' (inside a step) characteristics position in a
      * serial metadata buffer
      * <pre>
@@ -166,6 +171,22 @@ class VariableBase
      */
     void SetMemorySelection(const Box<Dims> &memorySelection);
 
+    /**
+     * Sets the requested accuracy for the next read operation.
+     * The actual accuracy after the read is provided in GetAccuracy()
+     */
+    void SetAccuracy(const adios2::Accuracy &a) noexcept;
+
+    /**
+     * Get the provided accuracy for the last read operation.
+     * Most operations provide data as it was written, meaning that
+     * error is reported as 0.0
+     */
+    adios2::Accuracy GetAccuracy() const noexcept;
+
+    /** Return the requested accuracy set by user with SetAccuracy */
+    adios2::Accuracy GetAccuracyRequested() const noexcept;
+
     size_t GetAvailableStepsStart() const;
 
     size_t GetAvailableStepsCount() const;
diff --git a/source/adios2/operator/OperatorFactory.cpp b/source/adios2/operator/OperatorFactory.cpp
index 3e9b0b419c..87139a92ff 100644
--- a/source/adios2/operator/OperatorFactory.cpp
+++ b/source/adios2/operator/OperatorFactory.cpp
@@ -31,6 +31,10 @@
 #include "adios2/operator/compress/CompressMGARDPlus.h"
 #endif
 
+#ifdef ADIOS2_HAVE_MGARD_MDR
+#include "adios2/operator/refactor/RefactorMDR.h"
+#endif
+
 #ifdef ADIOS2_HAVE_PNG
 #include "adios2/operator/compress/CompressPNG.h"
 #endif
@@ -74,6 +78,8 @@ std::string OperatorTypeToString(const Operator::OperatorType type)
         return "sz";
     case Operator::COMPRESS_ZFP:
         return "zfp";
+    case Operator::REFACTOR_MDR:
+        return "mdr";
     case Operator::PLUGIN_INTERFACE:
         return "plugin";
     default:
@@ -139,6 +145,12 @@ std::shared_ptr<Operator> MakeOperator(const std::string &type, const Params &pa
     {
 #ifdef ADIOS2_HAVE_ZFP
         ret = std::make_shared<compress::CompressZFP>(parameters);
+#endif
+    }
+    else if (typeLowerCase == "mdr")
+    {
+#ifdef ADIOS2_HAVE_MGARD_MDR
+        ret = std::make_shared<refactor::RefactorMDR>(parameters);
 #endif
     }
     else if (typeLowerCase == "plugin")
diff --git a/source/adios2/operator/OperatorFactory.h b/source/adios2/operator/OperatorFactory.h
index 50519f93a4..0439a81bd1 100644
--- a/source/adios2/operator/OperatorFactory.h
+++ b/source/adios2/operator/OperatorFactory.h
@@ -17,6 +17,8 @@ namespace adios2
 namespace core
 {
 
+std::string OperatorTypeToString(const Operator::OperatorType type);
+
 std::shared_ptr<Operator> MakeOperator(const std::string &type, const Params &parameters);
 
 size_t Decompress(const char *bufferIn, const size_t sizeIn, char *dataOut, MemorySpace memSpace,
diff --git a/source/adios2/operator/refactor/RefactorMDR.cpp b/source/adios2/operator/refactor/RefactorMDR.cpp
new file mode 100644
index 0000000000..b2a5518bfa
--- /dev/null
+++ b/source/adios2/operator/refactor/RefactorMDR.cpp
@@ -0,0 +1,536 @@
+/*
+ * Distributed under the OSI-approved Apache License, Version 2.0.  See
+ * accompanying file Copyright.txt for details.
+ *
+ * RefactorMDR.cpp :
+ *
+ *  Created on: Sep 14, 2023
+ *      Author: Norbert Podhorszki <pnorbert@ornl.gov>
+ */
+
+#include "RefactorMDR.h"
+#include "adios2/helper/adiosFunctions.h"
+#include "adios2/operator/compress/CompressNull.h"
+#include <assert.h>
+#include <cstring>
+
+#include <stdio.h>
+
+#include <mgard/compress_x.hpp>
+
+namespace adios2
+{
+namespace core
+{
+namespace refactor
+{
+
+RefactorMDR::RefactorMDR(const Params &parameters)
+: Operator("mdr", REFACTOR_MDR, "refactor", parameters)
+{
+    config.normalize_coordinates = false;
+    config.log_level = 1;
+    config.decomposition = mgard_x::decomposition_type::MultiDim;
+    config.domain_decomposition = mgard_x::domain_decomposition_type::MaxDim;
+    // config.domain_decomposition = mgard_x::domain_decomposition_type::Block;
+    // config.block_size = 64;
+
+    config.dev_type = mgard_x::device_type::AUTO;
+    config.prefetch = false;
+    // config.max_memory_footprint = max_memory_footprint;
+
+    config.lossless = mgard_x::lossless_type::Huffman_Zstd;
+    // mgard_x::lossless_type::Huffman, mgard_x::lossless_type::Huffman_Zstd,
+    // mgard_x::lossless_type::CPU_Lossless
+
+    // This should be changed to 4 for float later
+    config.total_num_bitplanes = 64;
+}
+
+size_t RefactorMDR::GetEstimatedSize(const size_t ElemCount, const size_t ElemSize,
+                                     const size_t ndims, const size_t *dims) const
+{
+    std::cout << "RefactorMDR::GetEstimatedSize() called \n";
+    mgard_x::data_type mtype =
+        (ElemSize == 8 ? mgard_x::data_type::Double : mgard_x::data_type::Float);
+    DataType datatype = (ElemSize == 8 ? DataType::Double : DataType::Float);
+    Dims dimsV(ndims);
+    for (size_t i = 0; i < ndims; ++i)
+    {
+        dimsV[i] = dims[i];
+    }
+    Dims convertedDims = ConvertDims(dimsV, datatype, 3);
+    mgard_x::DIM mgardDim = ndims;
+    std::vector<mgard_x::SIZE> mgardCount;
+    for (const auto &c : convertedDims)
+    {
+        mgardCount.push_back(c);
+    }
+
+    mgard_x::Config cfg(config); // copy of const config
+    if (mtype == mgard_x::data_type::Float)
+    {
+        cfg.total_num_bitplanes = 32;
+    }
+
+    auto s = mgard_x::MDR::MDRMaxOutputDataSize(mgardDim, mtype, mgardCount, config);
+
+    size_t sizeIn = helper::GetTotalSize(convertedDims, ElemSize);
+    std::cout << "RefactorMDR Estimated Max output size = " << s << " for input size = " << sizeIn
+              << std::endl;
+    return (size_t)s + 128; // count in the header
+};
+
+struct RefactorMDRHeader
+{
+    uint64_t ndims;
+    std::vector<uint64_t> dims; // ndims values
+    adios2::DataType type;
+    uint8_t mgard_version_major;
+    uint8_t mgard_version_minor;
+    uint8_t mgard_version_patch;
+    bool wasRefactored;
+    uint64_t metadataHeaderSize;
+    uint64_t metadataSize;
+    uint8_t nSubdomains;
+    uint8_t nLevels;
+    uint8_t nBitPlanes;
+};
+
+size_t RefactorMDR::Operate(const char *dataIn, const Dims &blockStart, const Dims &blockCount,
+                            const DataType type, char *bufferOut)
+{
+    const uint8_t bufferVersion = 1;
+    size_t bufferOutOffset = 0;
+
+    MakeCommonHeader(bufferOut, bufferOutOffset, bufferVersion);
+
+    Dims convertedDims = ConvertDims(blockCount, type, 3);
+
+    const size_t ndims = convertedDims.size();
+    if (ndims > 5)
+    {
+        helper::Throw<std::invalid_argument>("Operator", "RefactorMDR", "Operate",
+                                             "MGARD does not support data in " +
+                                                 std::to_string(ndims) + " dimensions");
+    }
+
+    // mgard V1 metadata
+    PutParameter(bufferOut, bufferOutOffset, ndims);
+    for (const auto &d : convertedDims)
+    {
+        PutParameter(bufferOut, bufferOutOffset, d);
+    }
+    PutParameter(bufferOut, bufferOutOffset, type);
+    PutParameter(bufferOut, bufferOutOffset, static_cast<uint8_t>(MGARD_VERSION_MAJOR));
+    PutParameter(bufferOut, bufferOutOffset, static_cast<uint8_t>(MGARD_VERSION_MINOR));
+    PutParameter(bufferOut, bufferOutOffset, static_cast<uint8_t>(MGARD_VERSION_PATCH));
+    // mgard V1 metadata end
+
+    // set type
+    mgard_x::data_type mgardType;
+    if (type == helper::GetDataType<float>())
+    {
+        mgardType = mgard_x::data_type::Float;
+        config.total_num_bitplanes = 32;
+    }
+    else if (type == helper::GetDataType<double>())
+    {
+        mgardType = mgard_x::data_type::Double;
+        config.total_num_bitplanes = 64;
+    }
+    else if (type == helper::GetDataType<std::complex<float>>())
+    {
+        mgardType = mgard_x::data_type::Float;
+        config.total_num_bitplanes = 32;
+    }
+    else if (type == helper::GetDataType<std::complex<double>>())
+    {
+        mgardType = mgard_x::data_type::Double;
+        config.total_num_bitplanes = 64;
+    }
+    else
+    {
+        helper::Throw<std::invalid_argument>("Operator", "RefactorMDR", "Operate",
+                                             "MGARD only supports float and double types");
+    }
+    // set type end
+
+    // set mgard style dim info
+    mgard_x::DIM mgardDim = ndims;
+    std::vector<mgard_x::SIZE> mgardCount;
+    for (const auto &c : convertedDims)
+    {
+        mgardCount.push_back(c);
+    }
+    // set mgard style dim info end
+
+    // input under this size will not be refactored
+    const size_t thresholdSize = 100000;
+
+    size_t sizeIn = helper::GetTotalSize(blockCount, helper::GetDataTypeSize(type));
+
+    if (sizeIn < thresholdSize)
+    {
+        /* disable compression and add marker in the header*/
+        PutParameter(bufferOut, bufferOutOffset, false);
+        headerSize = bufferOutOffset;
+        return 0;
+    }
+    PutParameter(bufferOut, bufferOutOffset, true);
+
+    mgard_x::MDR::RefactoredMetadata refactored_metadata;
+    mgard_x::MDR::RefactoredData refactored_data;
+    char *dataptr = const_cast<char *>(dataIn);
+    mgard_x::pin_memory(dataptr, sizeIn, config);
+    mgard_x::MDR::MDRefactor(mgardDim, mgardType, mgardCount, dataIn, refactored_metadata,
+                             refactored_data, config, false);
+    mgard_x::unpin_memory(dataptr, config);
+
+    size_t nbytes = SerializeRefactoredData(refactored_metadata, refactored_data,
+                                            bufferOut + bufferOutOffset, SIZE_MAX);
+
+    bufferOutOffset += nbytes;
+
+    return bufferOutOffset;
+}
+
+// return number of bytes written
+size_t RefactorMDR::SerializeRefactoredData(mgard_x::MDR::RefactoredMetadata &refactored_metadata,
+                                            mgard_x::MDR::RefactoredData &refactored_data,
+                                            char *buffer, size_t maxsize)
+{
+    size_t offset = 0;
+    size_t MDRHeaderSize = 0;
+
+    /* Metadata header */
+    const uint64_t metadata_header_size = refactored_metadata.header.size();
+    {
+        PutParameter<uint64_t>(buffer, offset, metadata_header_size);
+        std::memcpy(buffer + offset, refactored_metadata.header.data(), metadata_header_size);
+        offset += metadata_header_size;
+    }
+
+    /* Metadata */
+    std::vector<mgard_x::Byte> serialized_metadata = refactored_metadata.Serialize();
+    const uint64_t metadata_size = (uint64_t)serialized_metadata.size();
+    {
+        PutParameter<uint64_t>(buffer, offset, metadata_size);
+        std::memcpy(buffer + offset, serialized_metadata.data(), metadata_size);
+        offset += metadata_size;
+    }
+
+    std::cout << "MDR metadata seralized " << offset << " bytes. header = " << metadata_header_size
+              << " metadata = " << metadata_size << "\n";
+
+    /* 3D table of subdomain X level X bitplane offsets, not all of them has data */
+    uint8_t nSubdomains = refactored_metadata.metadata.size();
+    uint8_t nLevels = 0;
+    uint8_t nBitPlanes = 0;
+
+    for (size_t subdomain_id = 0; subdomain_id < refactored_metadata.metadata.size();
+         subdomain_id++)
+    {
+        nLevels = std::max(nLevels,
+                           (uint8_t)refactored_metadata.metadata[subdomain_id].level_sizes.size());
+        for (size_t level_idx = 0;
+             level_idx < refactored_metadata.metadata[subdomain_id].level_sizes.size(); level_idx++)
+        {
+            nBitPlanes = std::max(
+                nBitPlanes,
+                (uint8_t)refactored_metadata.metadata[subdomain_id].level_sizes[level_idx].size());
+        }
+    }
+
+    PutParameter<uint8_t>(buffer, offset, nSubdomains);
+    PutParameter<uint8_t>(buffer, offset, nLevels);
+    PutParameter<uint8_t>(buffer, offset, nBitPlanes);
+    uint64_t tableSize = (nSubdomains * nLevels * nBitPlanes);
+    uint64_t *table = (uint64_t *)(buffer + offset);
+    std::fill(table, table + tableSize, 0ULL);
+    offset += tableSize * sizeof(uint64_t);
+    MDRHeaderSize = offset;
+
+    /* Individual components of refactored data */
+    size_t nBlocks = 0;
+    uint64_t tableIdx = 0;
+    for (size_t subdomain_id = 0; subdomain_id < refactored_metadata.metadata.size();
+         subdomain_id++)
+    {
+        for (size_t level_idx = 0;
+             level_idx < refactored_metadata.metadata[subdomain_id].level_sizes.size(); level_idx++)
+        {
+            tableIdx = subdomain_id * nLevels * nBitPlanes + level_idx * nBitPlanes;
+            for (size_t bitplane_idx = 0;
+                 bitplane_idx <
+                 refactored_metadata.metadata[subdomain_id].level_sizes[level_idx].size();
+                 bitplane_idx++)
+            {
+                std::memcpy(buffer + offset,
+                            refactored_data.data[subdomain_id][level_idx][bitplane_idx],
+                            refactored_metadata.metadata[subdomain_id]
+                                .level_sizes[level_idx][bitplane_idx]);
+                table[tableIdx] = offset - MDRHeaderSize;
+                offset +=
+                    refactored_metadata.metadata[subdomain_id].level_sizes[level_idx][bitplane_idx];
+                ++nBlocks;
+                ++tableIdx;
+            }
+        }
+    }
+
+    std::cout << "MDR serialized " << offset << " bytes, MDR header size = " << MDRHeaderSize
+              << " subdomains = " << (size_t)nSubdomains << " levels = " << (size_t)nLevels
+              << " bitplanes = " << (size_t)nBitPlanes << " blocks = " << nBlocks << "\n";
+    return offset;
+}
+
+size_t RefactorMDR::GetHeaderSize() const { return headerSize; }
+
+size_t RefactorMDR::ReconstructV1(const char *bufferIn, const size_t sizeIn, char *dataOut)
+{
+    // Do NOT remove even if the buffer version is updated. Data might be still
+    // in lagacy formats. This function must be kept for backward compatibility.
+    // If a newer buffer format is implemented, create another function, e.g.
+    // ReconstructV1 and keep this function for reconstructing legacy data.
+
+    config.log_level = 1;
+    // double s = std::numeric_limits<double>::infinity();
+
+    size_t bufferInOffset = 0;
+
+    const size_t ndims = GetParameter<size_t, size_t>(bufferIn, bufferInOffset);
+    Dims blockCount(ndims);
+    for (size_t i = 0; i < ndims; ++i)
+    {
+        blockCount[i] = GetParameter<size_t, size_t>(bufferIn, bufferInOffset);
+    }
+    const DataType type = GetParameter<DataType>(bufferIn, bufferInOffset);
+    m_VersionInfo = " Data is compressed using MGARD Version " +
+                    std::to_string(GetParameter<uint8_t>(bufferIn, bufferInOffset)) + "." +
+                    std::to_string(GetParameter<uint8_t>(bufferIn, bufferInOffset)) + "." +
+                    std::to_string(GetParameter<uint8_t>(bufferIn, bufferInOffset)) +
+                    ". Please make sure a compatible version is used for decompression.";
+
+    const bool isRefactored = GetParameter<bool>(bufferIn, bufferInOffset);
+
+    if (!isRefactored)
+    {
+        // data was copied as is from this offset
+        headerSize += bufferInOffset;
+        m_AccuracyProvided = m_AccuracyRequested;
+        return 0;
+    }
+
+    size_t sizeOut = helper::GetTotalSize(blockCount, helper::GetDataTypeSize(type));
+
+    mgard_x::MDR::RefactoredMetadata refactored_metadata;
+
+    /* Metadata header */
+    {
+        const uint64_t metadata_header_size = GetParameter<uint64_t>(bufferIn, bufferInOffset);
+        refactored_metadata.header.resize(metadata_header_size);
+        std::memcpy(refactored_metadata.header.data(), bufferIn + bufferInOffset,
+                    metadata_header_size);
+        bufferInOffset += metadata_header_size;
+    }
+
+    /* Metadata */
+    {
+        const uint64_t metadata_size = GetParameter<uint64_t>(bufferIn, bufferInOffset);
+        std::vector<mgard_x::Byte> serialized_metadata;
+        serialized_metadata.resize(metadata_size);
+        std::memcpy(serialized_metadata.data(), bufferIn + bufferInOffset, metadata_size);
+        bufferInOffset += metadata_size;
+        refactored_metadata.Deserialize(serialized_metadata);
+        refactored_metadata.InitializeForReconstruction();
+    }
+
+    mgard_x::MDR::RefactoredData refactored_data;
+    refactored_data.InitializeForReconstruction(refactored_metadata);
+
+    const uint8_t nSubdomains = GetParameter<uint8_t>(bufferIn, bufferInOffset);
+    const uint8_t nLevels = GetParameter<uint8_t>(bufferIn, bufferInOffset);
+    const uint8_t nBitPlanes = GetParameter<uint8_t>(bufferIn, bufferInOffset);
+
+    const uint64_t tableSize = (nSubdomains * nLevels * nBitPlanes);
+    const uint64_t *table = (uint64_t *)(bufferIn + bufferInOffset);
+    bufferInOffset += tableSize * sizeof(uint64_t);
+
+    /*
+        Reconstruction
+    */
+    const char *componentData = bufferIn + bufferInOffset; // data pieces are here in memory
+
+    mgard_x::MDR::ReconstructedData reconstructed_data;
+    for (auto &metadata : refactored_metadata.metadata)
+    {
+        metadata.requested_tol = m_AccuracyRequested.error;
+        metadata.requested_s = m_AccuracyRequested.norm;
+    }
+    mgard_x::MDR::MDRequest(refactored_metadata, config);
+    /*for (auto &metadata : refactored_metadata.metadata)
+    {
+        metadata.PrintStatus();
+    }*/
+
+    bool first_reconstruction = true; // only will be needed with progressive reconstruction
+
+    // Assemble data pieces from buffer
+    {
+        uint64_t tableIdx;
+        int num_subdomains = refactored_metadata.metadata.size();
+        assert(nSubdomains == refactored_metadata.metadata.size());
+        for (int subdomain_id = 0; subdomain_id < num_subdomains; subdomain_id++)
+        {
+            mgard_x::MDR::MDRMetadata metadata = refactored_metadata.metadata[subdomain_id];
+            int num_levels = metadata.level_sizes.size();
+            for (int level_idx = 0; level_idx < num_levels; level_idx++)
+            {
+                assert(nLevels >= metadata.level_sizes.size());
+                tableIdx = subdomain_id * nLevels * nBitPlanes + level_idx * nBitPlanes;
+                int num_bitplanes = metadata.level_sizes[level_idx].size();
+                int loaded_bitplanes = metadata.loaded_level_num_bitplanes[level_idx];
+                int reqested_bitplanes = metadata.requested_level_num_bitplanes[level_idx];
+                assert(nBitPlanes >= metadata.requested_level_num_bitplanes[level_idx]);
+                for (int bitplane_idx = loaded_bitplanes; bitplane_idx < reqested_bitplanes;
+                     bitplane_idx++)
+                {
+
+                    uint64_t componentSize = refactored_metadata.metadata[subdomain_id]
+                                                 .level_sizes[level_idx][bitplane_idx];
+                    const mgard_x::Byte *cdata =
+                        reinterpret_cast<const mgard_x::Byte *>(componentData + table[tableIdx]);
+
+                    refactored_data.data[subdomain_id][level_idx][bitplane_idx] =
+                        const_cast<mgard_x::Byte *>(cdata);
+
+                    ++tableIdx;
+                }
+                if (first_reconstruction)
+                {
+                    // initialize level signs
+                    refactored_data.level_signs[subdomain_id][level_idx] =
+                        (bool *)malloc(sizeof(bool) * metadata.level_num_elems[level_idx]);
+                    std::memset(refactored_data.level_signs[subdomain_id][level_idx], 0,
+                                sizeof(bool) * metadata.level_num_elems[level_idx]);
+                }
+            }
+        }
+    }
+
+    /* Initialize reconstructed data manually here to force using
+        user allocated memory for the final result
+    */
+    if (false)
+    {
+        reconstructed_data.Initialize(1);
+        reconstructed_data.data[0] = reinterpret_cast<mgard_x::Byte *>(dataOut);
+        std::memset(reconstructed_data.data[0], 0, sizeOut);
+        std::vector<mgard_x::SIZE> offsets = std::vector<mgard_x::SIZE>(ndims, 0);
+        std::vector<mgard_x::SIZE> shape = std::vector<mgard_x::SIZE>();
+        for (const auto &c : blockCount)
+        {
+            shape.push_back(static_cast<mgard_x::SIZE>(c));
+        }
+        reconstructed_data.offset[0] = offsets;
+        reconstructed_data.shape[0] = shape;
+    }
+
+    mgard_x::MDR::MDReconstruct(refactored_metadata, refactored_data, reconstructed_data, config,
+                                false);
+
+    size_t refactoredSize =
+        helper::GetTotalSize(reconstructed_data.shape[0], helper::GetDataTypeSize(type));
+    assert(sizeOut == refactoredSize);
+    std::memcpy(dataOut, reconstructed_data.data[0], sizeOut);
+
+    first_reconstruction = false;
+
+    /*for (int subdomain_id = 0; subdomain_id < reconstructed_data.data.size(); subdomain_id++)
+    {
+        std::cout << "reconstructed_data " << subdomain_id << " : offset(";
+        for (auto n : reconstructed_data.offset[subdomain_id])
+        {
+            std::cout << n << " ";
+        }
+        std::cout << ") shape(";
+        for (auto n : reconstructed_data.shape[subdomain_id])
+        {
+            std::cout << n << " ";
+        }
+        std::cout << ")\n";
+    }*/
+
+    if (type == DataType::FloatComplex || type == DataType::Float)
+    {
+        m_AccuracyProvided.error =
+            std::max(m_AccuracyRequested.error, adios2::ops::mdr::FLOAT_ROUNDING_ERROR_LIMIT);
+    }
+    else
+    {
+        m_AccuracyProvided.error =
+            std::max(m_AccuracyRequested.error, adios2::ops::mdr::DOUBLE_ROUNDING_ERROR_LIMIT);
+    }
+    m_AccuracyProvided.norm = m_AccuracyRequested.norm;
+    m_AccuracyProvided.relative = false; // should be m_AccuracyRequested.relative
+
+    if (type == DataType::FloatComplex || type == DataType::DoubleComplex)
+    {
+        sizeOut /= 2;
+    }
+    return sizeOut;
+}
+
+size_t RefactorMDR::InverseOperate(const char *bufferIn, const size_t sizeIn, char *dataOut)
+{
+
+    for (auto &p : m_Parameters)
+    {
+        std::cout << "User parameter " << p.first << " = " << p.second << std::endl;
+        const std::string key = helper::LowerCase(p.first);
+        // const std::string value = helper::LowerCase(p.second);
+        if (key == "accuracy")
+        {
+            m_AccuracyRequested.error =
+                helper::StringTo<double>(p.second, " in Parameter key=" + key);
+            std::cout << "Accuracy error set from Parameter to " << m_AccuracyRequested.error;
+        }
+    }
+
+    size_t bufferInOffset = 1; // skip operator type
+    const uint8_t bufferVersion = GetParameter<uint8_t>(bufferIn, bufferInOffset);
+    bufferInOffset += 2; // skip two reserved bytes
+    headerSize = bufferInOffset;
+
+    if (bufferVersion == 1)
+    {
+        return ReconstructV1(bufferIn + bufferInOffset, sizeIn - bufferInOffset, dataOut);
+    }
+    /*else if (bufferVersion == 2)
+    {
+        // TODO: if a Version 2 mgard buffer is being implemented, put it here
+        // and keep the ReconstructV1 routine for backward compatibility
+    }*/
+    else
+    {
+        helper::Throw<std::runtime_error>("Operator", "RefactorMDR", "InverseOperate",
+                                          "invalid mgard buffer version" + bufferVersion);
+    }
+
+    return 0;
+}
+
+bool RefactorMDR::IsDataTypeValid(const DataType type) const
+{
+    if (type == DataType::Double || type == DataType::Float || type == DataType::DoubleComplex ||
+        type == DataType::FloatComplex)
+    {
+        return true;
+    }
+    return false;
+}
+
+} // end namespace compress
+} // end namespace core
+} // end namespace adios2
diff --git a/source/adios2/operator/refactor/RefactorMDR.h b/source/adios2/operator/refactor/RefactorMDR.h
new file mode 100644
index 0000000000..99a39a82a2
--- /dev/null
+++ b/source/adios2/operator/refactor/RefactorMDR.h
@@ -0,0 +1,85 @@
+/*
+ * Distributed under the OSI-approved Apache License, Version 2.0.  See
+ * accompanying file Copyright.txt for details.
+ *
+ * RefactorMDR.h :
+ *
+ *  Created on: Sep 14, 2023
+ *      Author: Norbert Podhorszki <pnorbert@ornl.gov>
+ */
+
+#ifndef ADIOS2_OPERATOR_COMPRESS_REFACTORMDR_H_
+#define ADIOS2_OPERATOR_COMPRESS_REFACTORMDR_H_
+
+#include "adios2/core/Operator.h"
+
+#include <mgard/MGARDConfig.hpp>
+#include <mgard/mdr_x.hpp>
+
+namespace adios2
+{
+namespace core
+{
+namespace refactor
+{
+
+class RefactorMDR : public Operator
+{
+
+public:
+    RefactorMDR(const Params &parameters);
+
+    ~RefactorMDR() = default;
+
+    /**
+     * @param dataIn
+     * @param blockStart
+     * @param blockCount
+     * @param type
+     * @param bufferOut
+     * @return size of compressed buffer
+     */
+    size_t Operate(const char *dataIn, const Dims &blockStart, const Dims &blockCount,
+                   const DataType type, char *bufferOut) final;
+
+    /**
+     * @param bufferIn
+     * @param sizeIn
+     * @param dataOut
+     * @return size of decompressed buffer
+     */
+    size_t InverseOperate(const char *bufferIn, const size_t sizeIn, char *dataOut) final;
+
+    bool IsDataTypeValid(const DataType type) const final;
+
+    size_t GetHeaderSize() const;
+
+    size_t GetEstimatedSize(const size_t ElemCount, const size_t ElemSize, const size_t ndims,
+                            const size_t *dims) const;
+
+private:
+    size_t headerSize = 0;
+    /**
+     * Decompress function for V1 buffer. Do NOT remove even if the buffer
+     * version is updated. Data might be still in lagacy formats. This function
+     * must be kept for backward compatibility
+     * @param bufferIn : compressed data buffer (V1 only)
+     * @param sizeIn : number of bytes in bufferIn
+     * @param dataOut : decompressed data buffer
+     * @return : number of bytes in dataOut
+     */
+    size_t ReconstructV1(const char *bufferIn, const size_t sizeIn, char *dataOut);
+    size_t SerializeRefactoredData(mgard_x::MDR::RefactoredMetadata &refactored_metadata,
+                                   mgard_x::MDR::RefactoredData &refactored_data, char *buffer,
+                                   size_t maxsize);
+
+    std::string m_VersionInfo;
+
+    mgard_x::Config config;
+};
+
+} // end namespace compress
+} // end namespace core
+} // end namespace adios2
+
+#endif /* ADIOS2_OPERATOR_COMPRESS_REFACTORMDR_H_ */
diff --git a/source/adios2/toolkit/format/bp5/BP5Deserializer.cpp b/source/adios2/toolkit/format/bp5/BP5Deserializer.cpp
index 6e6ea7e806..9e5cf4067b 100644
--- a/source/adios2/toolkit/format/bp5/BP5Deserializer.cpp
+++ b/source/adios2/toolkit/format/bp5/BP5Deserializer.cpp
@@ -1778,12 +1778,29 @@ void BP5Deserializer::FinalizeGet(const ReadRequest &Read, const bool freeAddr)
             DestSize *= writer_meta_base->Count[dim + Read.BlockID * writer_meta_base->Dims];
         }
         decompressBuffer.resize(DestSize);
+
+        // Get the operator of the variable if exists or create one
+        std::shared_ptr<Operator> op = nullptr;
+        VariableBase *VB = static_cast<VariableBase *>(((struct BP5VarRec *)Req.VarRec)->Variable);
+        if (!VB->m_Operations.empty())
+        {
+            op = VB->m_Operations[0];
+        }
+        else
+        {
+            Operator::OperatorType compressorType =
+                static_cast<Operator::OperatorType>(IncomingData[0]);
+            op = MakeOperator(OperatorTypeToString(compressorType), {});
+        }
+        op->SetAccuracy(VB->GetAccuracyRequested());
+
         {
             std::lock_guard<std::mutex> lockGuard(mutexDecompress);
             core::Decompress(
                 IncomingData,
                 ((MetaArrayRecOperator *)writer_meta_base)->DataBlockSize[Read.BlockID],
-                decompressBuffer.data(), Req.MemSpace);
+                decompressBuffer.data(), Req.MemSpace, op);
+            VB->m_AccuracyProvided = op->GetAccuracy();
         }
         IncomingData = decompressBuffer.data();
         VirtualIncomingData = IncomingData;
diff --git a/source/adios2/toolkit/format/bp5/BP5Serializer.cpp b/source/adios2/toolkit/format/bp5/BP5Serializer.cpp
index d7145bd7fb..84fcd33fc3 100644
--- a/source/adios2/toolkit/format/bp5/BP5Serializer.cpp
+++ b/source/adios2/toolkit/format/bp5/BP5Serializer.cpp
@@ -837,7 +837,8 @@ void BP5Serializer::Marshal(void *Variable, const char *Name, const DataType Typ
                 if (Offsets)
                     tmpOffsets.push_back(Offsets[i]);
             }
-            size_t AllocSize = ElemCount * ElemSize + 100;
+            size_t AllocSize =
+                VB->m_Operations[0]->GetEstimatedSize(ElemCount, ElemSize, DimCount, Count);
             BufferV::BufferPos pos = CurDataBuffer->Allocate(AllocSize, ElemSize);
             char *CompressedData = (char *)GetPtr(pos.bufferIdx, pos.posInBuffer);
             DataOffset = m_PriorDataBufferSizeTotal + pos.globalPos;
diff --git a/source/utils/bpls/bpls.cpp b/source/utils/bpls/bpls.cpp
index c696690fa0..7b80fd1c34 100644
--- a/source/utils/bpls/bpls.cpp
+++ b/source/utils/bpls/bpls.cpp
@@ -83,6 +83,7 @@ std::string format;           // format string for one data element (e.g. %6.2f)
 std::string transport_params; // Transport parameters (e.g. "Library=stdio,verbose=3")
 std::string engine_name;      // Engine name (e.g. "BP5")
 std::string engine_params;    // Engine parameters (e.g. "SelectSteps=0:5:2")
+std::string accuracy_def;     // Accuracy definition (e.g. "accuracy="0.0,0.0,rel")
 
 // Flags from arguments or defaults
 bool dump; // dump data not just list info(flag == 1)
@@ -102,6 +103,8 @@ bool hidden_attrs;       // show hidden attrs in BP file
 int hidden_attrs_flag;   // to be passed on in option struct
 bool show_decomp;        // show decomposition of arrays
 bool show_version;       // print binary version info of file before work
+adios2::Accuracy accuracy;
+bool accuracyWasSet = false;
 
 // other global variables
 char *prgname; /* argv[0] */
@@ -169,6 +172,10 @@ void display_help()
            "file\n"
            "  --decomp    | -D           Show decomposition of variables as layed "
            "out in file\n"
+           "  --error     | -X string    Specify read accuracy (error,norm,rel|abs)\n"
+           "                             e.g. error=\"0.0,0.0,abs\"\n"
+           "                             L2 norm = 0.0, Linf = inf\n"
+
            "  --transport-parameters | -T         Specify File transport "
            "parameters\n"
            "                                      e.g. \"Library=stdio\"\n"
@@ -636,6 +643,9 @@ int bplsMain(int argc, char *argv[])
                            "Print version information (add -verbose for additional"
                            " information)");
     arg.AddBooleanArgument("-V", &show_version, "");
+    arg.AddArgument("--error", argT::SPACE_ARGUMENT, &accuracy_def,
+                    "| -X string    Specify read accuracy (error,norm,rel|abs)");
+    arg.AddArgument("-X", argT::SPACE_ARGUMENT, &accuracy_def, "");
     arg.AddArgument("--transport-parameters", argT::SPACE_ARGUMENT, &transport_params,
                     "| -T string    Specify File transport parameters manually");
     arg.AddArgument("-T", argT::SPACE_ARGUMENT, &transport_params, "");
@@ -705,6 +715,10 @@ int bplsMain(int argc, char *argv[])
     if (attrsonly)
         listattrs = true;
 
+    retval = parseAccuracy();
+    if (retval)
+        return retval;
+
     if (verbose > 1)
         printSettings();
 
@@ -1273,6 +1287,7 @@ int printVariableInfo(core::Engine *fp, core::IO *io, core::Variable<T> *variabl
 
     if (dump && !show_decomp)
     {
+        variable->SetAccuracy(accuracy);
         // print variable content
         if (variable->m_ShapeID == ShapeID::LocalArray)
         {
@@ -2210,6 +2225,14 @@ int readVar(core::Engine *fp, core::IO *io, core::Variable<T> *variable)
         }
     } // end while sum < nelems
     print_endline();
+
+    if (accuracyWasSet)
+    {
+        adios2::Accuracy a = variable->GetAccuracy();
+        std::cout << "Read accuracy was (error = " << a.error << ", norm = " << a.norm << ", "
+                  << (a.relative ? "rel" : "abs") << ")\n";
+    }
+
     return 0;
 }
 
@@ -3802,6 +3825,61 @@ void print_decomp_singlestep(core::Engine *fp, core::IO *io, core::Variable<T> *
     }
 }
 
+int parseAccuracy()
+{
+    if (accuracy_def.empty())
+        return 0;
+
+    // Vector of string to save tokens
+    std::vector<std::string> tokens;
+
+    // stringstream class check1
+    std::stringstream ss(accuracy_def);
+
+    std::string intermediate;
+
+    // Tokenizing w.r.t. space ','
+    while (getline(ss, intermediate, ','))
+    {
+        tokens.push_back(intermediate);
+    }
+
+    if (tokens.size() != 3)
+    {
+        std::cout
+            << "ERROR: --error definition needs 3 values: error, norm, and relative|absolute\n"
+            << " error >= 0.0, norm >= 0.0 or inf, third arg is \"rel\" or \"abs\"\n";
+        return 1;
+    }
+
+    accuracy.error = adios2::helper::StringTo<double>(tokens[0], "error value for accuracy");
+    std::string normval = adios2::helper::LowerCase(tokens[2]);
+    if (normval == "inf")
+        accuracy.norm = std::numeric_limits<double>::infinity();
+    else
+        accuracy.norm = adios2::helper::StringTo<double>(tokens[1], "norm value for accuracy");
+    std::string relval = adios2::helper::LowerCase(tokens[2]);
+    if (relval == "rel")
+        accuracy.relative = true;
+    else if (relval == "abs")
+        accuracy.relative = false;
+    else
+    {
+        std::cout << "ERROR: --error third value must be \"rel\" or \"abs\" but it was \""
+                  << tokens[2] << "\"\n";
+        return 1;
+    }
+
+    accuracyWasSet = true;
+    if (verbose > 0)
+    {
+        std::cout << "Read accuracy is set to (error = " << accuracy.error
+                  << ", norm = " << accuracy.norm << ", " << (accuracy.relative ? "rel" : "abs")
+                  << ")\n";
+    }
+    return 0;
+}
+
 // parse a string "0, 3; 027" into an integer array
 // of [0,3,27]
 // exits if parsing failes
diff --git a/source/utils/bpls/bpls.h b/source/utils/bpls/bpls.h
index f0a715a94f..683c582b4d 100644
--- a/source/utils/bpls/bpls.h
+++ b/source/utils/bpls/bpls.h
@@ -57,6 +57,7 @@ char *mystrndup(const char *s, size_t n);
 void init_globals();
 void processDimSpecs();
 void parseDimSpec(const std::string &str, int64_t *dims);
+int parseAccuracy();
 int compile_regexp_masks(void);
 void printSettings(void);
 int doList(const char *path);
diff --git a/testing/adios2/engine/bp/CMakeLists.txt b/testing/adios2/engine/bp/CMakeLists.txt
index 03479c7edc..8c5b90606b 100644
--- a/testing/adios2/engine/bp/CMakeLists.txt
+++ b/testing/adios2/engine/bp/CMakeLists.txt
@@ -184,6 +184,11 @@ gtest_add_tests_helper(ReadMultithreaded MPI_NONE BP Engine.BP. .BP5
   WORKING_DIRECTORY ${BP5_DIR} EXTRA_ARGS "BP5"
 )
 
+# Only a single test is enough, pick the latest engine
+gtest_add_tests_helper(AccuracyDefaults MPI_NONE BP Engine.BP. .BP5
+  WORKING_DIRECTORY ${BP5_DIR} EXTRA_ARGS "BP5"
+)
+
 # BP3 only for now
 gtest_add_tests_helper(WriteNull MPI_ALLOW BP Engine.BP. .BP3
   WORKING_DIRECTORY ${BP3_DIR} EXTRA_ARGS "BP3"
@@ -191,6 +196,7 @@ gtest_add_tests_helper(WriteNull MPI_ALLOW BP Engine.BP. .BP3
 
 # BP4 and BP5 but NOT BP3
 bp4_bp5_gtest_add_tests_helper(WriteAppendReadADIOS2 MPI_ALLOW)
+bp4_bp5_gtest_add_tests_helper(JoinedArray MPI_ALLOW)
 
 # BP4 only for now
 # gtest_add_tests_helper(WriteAppendReadADIOS2 MPI_ALLOW BP Engine.BP. .BP4
@@ -202,9 +208,9 @@ gtest_add_tests_helper(StepsInSituGlobalArray MPI_ALLOW BP Engine.BP. .BP4
 gtest_add_tests_helper(StepsInSituLocalArray MPI_ALLOW BP Engine.BP. .BP4
   WORKING_DIRECTORY ${BP4_DIR} EXTRA_ARGS "BP4"
 )
-gtest_add_tests_helper(JoinedArray MPI_ALLOW BP Engine.BP. .BP4
-  WORKING_DIRECTORY ${BP4_DIR} EXTRA_ARGS "BP4"
-)
+#gtest_add_tests_helper(JoinedArray MPI_ALLOW BP Engine.BP. .BP4
+#  WORKING_DIRECTORY ${BP4_DIR} EXTRA_ARGS "BP4"
+#)
 
 # InquireVaribleException only for BP4 because BP5 works differently
 gtest_add_tests_helper(InquireVariableException MPI_ALLOW BP Engine.BP. .BP4
diff --git a/testing/adios2/engine/bp/TestBPAccuracyDefaults.cpp b/testing/adios2/engine/bp/TestBPAccuracyDefaults.cpp
new file mode 100644
index 0000000000..d0493cdc22
--- /dev/null
+++ b/testing/adios2/engine/bp/TestBPAccuracyDefaults.cpp
@@ -0,0 +1,153 @@
+/*
+ * Distributed under the OSI-approved Apache License, Version 2.0.  See
+ * accompanying file Copyright.txt for details.
+ *
+ * TestBPAccuracyDefaults.cpp :
+ *
+ *  Created on: Dec 5, 2023
+ *      Author: Norbert Podhorszki
+ */
+
+#include <cstdint>
+#include <cstring>
+
+#include <iostream>
+#include <numeric> //std::iota
+#include <stdexcept>
+
+#include <adios2.h>
+
+#include <gtest/gtest.h>
+
+std::string engineName;       // comes from command line
+std::string engineParameters; // comes from command line
+
+class AccuracyTests : public ::testing::Test
+{
+public:
+    AccuracyTests() = default;
+};
+
+// Check if SetAccuracy/GetAccuracy default behavior works
+TEST_F(AccuracyTests, DefaultAccuracy)
+{
+    const std::string fname("DefaultAccuracy.bp");
+
+    int mpiRank = 0, mpiSize = 1;
+
+    const std::size_t Nx = 10;
+
+#if ADIOS2_USE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
+    MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
+#endif
+
+    std::vector<int64_t> localData(Nx);
+    std::iota(localData.begin(), localData.end(), mpiRank * Nx);
+
+#if ADIOS2_USE_MPI
+    adios2::ADIOS adios(MPI_COMM_WORLD);
+#else
+    adios2::ADIOS adios;
+#endif
+    {
+        adios2::IO io = adios.DeclareIO("DefaultAccuracyWrite");
+        if (!engineName.empty())
+        {
+            io.SetEngine(engineName);
+        }
+        if (!engineParameters.empty())
+        {
+            io.SetParameters(engineParameters);
+        }
+
+        adios2::Variable<int64_t> var =
+            io.DefineVariable<int64_t>("range", {static_cast<std::size_t>(Nx * mpiSize)},
+                                       {static_cast<std::size_t>(Nx * mpiRank)}, {Nx});
+
+        adios2::Accuracy accuracyRequested = {0.001, std::numeric_limits<double>::infinity(), true};
+        var.SetAccuracy(accuracyRequested);
+        adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write);
+
+        bpWriter.Put<int64_t>("range", localData.data());
+        bpWriter.Close();
+
+        auto accuracyGot = var.GetAccuracy();
+        assert(accuracyGot.error == 0.0); // no error whatsoever
+        assert(accuracyGot.norm == accuracyRequested.norm);
+        assert(accuracyGot.relative == accuracyRequested.relative);
+    }
+    // Reader
+    {
+        adios2::IO io = adios.DeclareIO("DefaultAccuracyRead");
+        if (!engineName.empty())
+        {
+            io.SetEngine(engineName);
+        }
+        if (!engineParameters.empty())
+        {
+            io.SetParameters(engineParameters);
+        }
+
+        adios2::Engine bpReader = io.Open(fname, adios2::Mode::ReadRandomAccess);
+        adios2::Variable<int64_t> varRange = io.InquireVariable<int64_t>("range");
+        adios2::Accuracy accuracyRequested = {0.001, 1.17, false};
+        varRange.SetAccuracy(accuracyRequested);
+
+        const std::size_t gNx = static_cast<std::size_t>(Nx * mpiSize);
+        std::vector<int64_t> globalData(gNx);
+        bpReader.Get(varRange, globalData);
+        bpReader.PerformGets();
+
+        auto accuracyGot = varRange.GetAccuracy();
+        assert(accuracyGot.error == 0.0); // no error whatsoever
+        assert(accuracyGot.norm == accuracyRequested.norm);
+        assert(accuracyGot.relative == accuracyRequested.relative);
+
+        std::vector<int64_t> iStartEndData;
+        iStartEndData.reserve(gNx); // maximum possible
+
+        for (size_t i = 1; i < gNx; ++i)
+        {
+            varRange.SetSelection({{i}, {gNx - i}});
+
+            bpReader.Get("range", iStartEndData);
+            bpReader.PerformGets();
+
+            for (size_t j = i; j < gNx; ++j)
+            {
+                EXPECT_EQ(globalData[j], iStartEndData[j - i]);
+            }
+        }
+        bpReader.Close();
+    }
+}
+
+int main(int argc, char **argv)
+{
+#if ADIOS2_USE_MPI
+    int provided;
+
+    // MPI_THREAD_MULTIPLE is only required if you enable the SST MPI_DP
+    MPI_Init_thread(nullptr, nullptr, MPI_THREAD_MULTIPLE, &provided);
+#endif
+
+    int result;
+    ::testing::InitGoogleTest(&argc, argv);
+
+    if (argc > 1)
+    {
+        engineName = std::string(argv[1]);
+    }
+    if (argc > 2)
+    {
+        engineParameters = std::string(argv[2]);
+    }
+    result = RUN_ALL_TESTS();
+
+#if ADIOS2_USE_MPI
+    MPI_Finalize();
+#endif
+
+    return result;
+}
diff --git a/testing/adios2/engine/bp/operations/CMakeLists.txt b/testing/adios2/engine/bp/operations/CMakeLists.txt
index 61a4544b01..b8c2a04362 100644
--- a/testing/adios2/engine/bp/operations/CMakeLists.txt
+++ b/testing/adios2/engine/bp/operations/CMakeLists.txt
@@ -60,6 +60,12 @@ if(ADIOS2_HAVE_MGARD)
       target_link_libraries(${tgt} CUDA::cudart)
     endforeach()
   endif()
+
+  if(ADIOS2_HAVE_MGARD_MDR)
+    gtest_add_tests_helper(WriteReadMGARDMDR MPI_ALLOW BP Engine.BP. .BP5
+    WORKING_DIRECTORY ${BP5_DIR} EXTRA_ARGS "BP5"
+    )
+  endif()
 endif()
 
 if(ADIOS2_HAVE_BZip2)
diff --git a/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDMDR.cpp b/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDMDR.cpp
new file mode 100644
index 0000000000..084e21631b
--- /dev/null
+++ b/testing/adios2/engine/bp/operations/TestBPWriteReadMGARDMDR.cpp
@@ -0,0 +1,227 @@
+/*
+ * Distributed under the OSI-approved Apache License, Version 2.0.  See
+ * accompanying file Copyright.txt for details.
+ */
+#include <cstdint>
+#include <cstring>
+
+#include <algorithm>
+#include <iostream> //std::cout
+#include <numeric>  //std::iota
+#include <stdexcept>
+
+#include <adios2.h>
+
+#include <gtest/gtest.h>
+
+std::string engineName; // comes from command line
+
+class BPWriteReadMGARDMDR : public ::testing::TestWithParam<std::string>
+{
+public:
+    BPWriteReadMGARDMDR() = default;
+    virtual void SetUp(){};
+    virtual void TearDown(){};
+};
+
+TEST_F(BPWriteReadMGARDMDR, BPWRMGARD1D)
+{
+    // Refactor a dataset with MDR, then
+    // read back with various accuracies
+    const std::string fname("BPWRMGARDMDR1D.bp");
+
+    int mpiRank = 0, mpiSize = 1;
+    const size_t Nx = 30000; // 100k minimum data size for MDR
+    const size_t NSteps = 1;
+
+    std::vector<float> r32s(Nx);
+    std::vector<double> r64s(Nx);
+
+    const double pi = 3.141592653589793238462643383279;
+    const double twopi = 2 * pi;
+    const double value = (2 * pi) / double(Nx);
+    for (size_t x = 0; x < Nx; ++x)
+    {
+        r64s[x] = std::cos(value * x);
+        r32s[x] = (float)r64s[x];
+    }
+
+#if ADIOS2_USE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
+    MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
+#endif
+
+#if ADIOS2_USE_MPI
+    adios2::ADIOS adios(MPI_COMM_WORLD);
+#else
+    adios2::ADIOS adios;
+#endif
+    {
+        adios2::IO io = adios.DeclareIO("WriteIO");
+
+        if (!engineName.empty())
+        {
+            io.SetEngine(engineName);
+        }
+
+        const adios2::Dims shape{static_cast<size_t>(Nx * mpiSize)};
+        const adios2::Dims start{static_cast<size_t>(Nx * mpiRank)};
+        const adios2::Dims count{Nx};
+
+        auto var_r32 = io.DefineVariable<float>("r32", shape, start, count, adios2::ConstantDims);
+        auto var_r64 = io.DefineVariable<double>("r64", shape, start, count, adios2::ConstantDims);
+        auto var_raw64 =
+            io.DefineVariable<double>("raw64", shape, start, count, adios2::ConstantDims);
+
+        // add operations
+        adios2::Operator mgardOp = adios.DefineOperator("mgardCompressor", adios2::ops::MDR);
+
+        var_r32.AddOperation(mgardOp, {});
+        var_r64.AddOperation(mgardOp, {});
+
+        adios2::Engine bpWriter = io.Open(fname, adios2::Mode::Write);
+
+        for (size_t step = 0; step < NSteps; ++step)
+        {
+            bpWriter.BeginStep();
+            bpWriter.Put<float>("r32", r32s.data());
+            bpWriter.Put<double>("r64", r64s.data());
+            bpWriter.Put<double>("raw64", r64s.data());
+            bpWriter.EndStep();
+        }
+
+        bpWriter.Close();
+    }
+
+    {
+        adios2::IO io = adios.DeclareIO("ReadIO");
+
+        if (!engineName.empty())
+        {
+            io.SetEngine(engineName);
+        }
+
+        adios2::Engine bpReader = io.Open(fname, adios2::Mode::Read);
+        size_t t = 0;
+        while (bpReader.BeginStep() == adios2::StepStatus::OK)
+        {
+            auto var_r32 = io.InquireVariable<float>("r32");
+            EXPECT_TRUE(var_r32);
+            ASSERT_EQ(var_r32.ShapeID(), adios2::ShapeID::GlobalArray);
+            ASSERT_EQ(var_r32.Steps(), NSteps);
+            ASSERT_EQ(var_r32.Shape()[0], mpiSize * Nx);
+
+            auto var_r64 = io.InquireVariable<double>("r64");
+            EXPECT_TRUE(var_r64);
+            ASSERT_EQ(var_r64.ShapeID(), adios2::ShapeID::GlobalArray);
+            ASSERT_EQ(var_r64.Steps(), NSteps);
+            ASSERT_EQ(var_r64.Shape()[0], mpiSize * Nx);
+
+            const adios2::Dims start{mpiRank * Nx};
+            const adios2::Dims count{Nx};
+            const adios2::Box<adios2::Dims> sel(start, count);
+            var_r32.SetSelection(sel);
+            var_r64.SetSelection(sel);
+
+            std::vector<double> errors = {0.0, 0.0000001, 0.0001, 0.1};
+            for (auto error : errors)
+            {
+                std::cout << "===== Read with error = " << error << std::endl;
+
+                std::vector<float> read32s;
+                std::vector<double> read64s;
+
+                adios2::Accuracy accuracyRequested = {
+                    error, std::numeric_limits<double>::infinity(), false};
+                var_r32.SetAccuracy(accuracyRequested);
+                var_r64.SetAccuracy(accuracyRequested);
+                bpReader.Get(var_r32, read32s, adios2::Mode::Sync);
+                bpReader.Get(var_r64, read64s, adios2::Mode::Sync);
+
+                auto accuracyGot32 = var_r32.GetAccuracy();
+                assert(accuracyGot32.error <=
+                       std::max(accuracyRequested.error,
+                                adios2::ops::mdr::FLOAT_ROUNDING_ERROR_LIMIT));
+                assert(accuracyGot32.norm == accuracyRequested.norm);
+                assert(accuracyGot32.relative == accuracyRequested.relative);
+
+                auto accuracyGot64 = var_r64.GetAccuracy();
+                assert(accuracyGot64.error <=
+                       std::max(accuracyRequested.error,
+                                adios2::ops::mdr::DOUBLE_ROUNDING_ERROR_LIMIT));
+                assert(accuracyGot64.norm == accuracyRequested.norm);
+                assert(accuracyGot64.relative == accuracyRequested.relative);
+
+                double maxDiff = 0.0, relativeMaxDiff = 0.0;
+                size_t maxDiffPos = 0;
+
+                for (size_t i = 0; i < Nx; ++i)
+                {
+                    double diff = std::abs(r64s[i] - read64s[i]);
+                    if (diff > maxDiff)
+                    {
+                        maxDiff = diff;
+                        maxDiffPos = i;
+                    }
+                }
+
+                auto r64s_Max = std::max_element(r64s.begin(), r64s.end());
+                relativeMaxDiff = maxDiff / *r64s_Max;
+                std::cout << "double array: Relative Max Diff " << relativeMaxDiff << " Max Diff "
+                          << maxDiff << " requested error " << error << " result error "
+                          << accuracyGot64.error << " max diff pos " << maxDiffPos << " orig value "
+                          << r64s[maxDiffPos] << " read value " << read64s[maxDiffPos] << "\n";
+                ASSERT_LE(maxDiff, accuracyGot64.error);
+
+                for (size_t i = 0; i < Nx; ++i)
+                {
+                    double diff = std::abs(r32s[i] - read32s[i]);
+                    if (diff > maxDiff)
+                    {
+                        maxDiff = diff;
+                        maxDiffPos = i;
+                    }
+                }
+
+                auto r32s_Max = std::max_element(r32s.begin(), r32s.end());
+                relativeMaxDiff = maxDiff / *r32s_Max;
+
+                std::cout << "float array: Relative Max Diff " << relativeMaxDiff << " Max Diff "
+                          << maxDiff << " requested error " << error << " result error "
+                          << accuracyGot32.error << " max diff pos " << maxDiffPos << " orig value "
+                          << r32s[maxDiffPos] << " read value " << read32s[maxDiffPos] << "\n";
+                ASSERT_LE(maxDiff, accuracyGot32.error);
+            }
+            bpReader.EndStep();
+            ++t;
+        }
+
+        EXPECT_EQ(t, NSteps);
+
+        bpReader.Close();
+    }
+}
+
+int main(int argc, char **argv)
+{
+#if ADIOS2_USE_MPI
+    int provided;
+
+    // MPI_THREAD_MULTIPLE is only required if you enable the SST MPI_DP
+    MPI_Init_thread(nullptr, nullptr, MPI_THREAD_MULTIPLE, &provided);
+#endif
+
+    int result;
+    ::testing::InitGoogleTest(&argc, argv);
+    if (argc > 1)
+    {
+        engineName = std::string(argv[1]);
+    }
+    result = RUN_ALL_TESTS();
+
+#if ADIOS2_USE_MPI
+    MPI_Finalize();
+#endif
+
+    return result;
+}
diff --git a/testing/utils/cwriter/TestUtilsCWriter.bplsh.expected.txt b/testing/utils/cwriter/TestUtilsCWriter.bplsh.expected.txt
index d9461a11fa..24bd5cd4c4 100644
--- a/testing/utils/cwriter/TestUtilsCWriter.bplsh.expected.txt
+++ b/testing/utils/cwriter/TestUtilsCWriter.bplsh.expected.txt
@@ -28,6 +28,9 @@ The time dimension is the first dimension then.
                                instead of the default. E.g. "%6.3f"
   --hidden_attrs             Show hidden ADIOS attributes in the file
   --decomp    | -D           Show decomposition of variables as layed out in file
+  --error     | -X string    Specify read accuracy (error,norm,rel|abs)
+                             e.g. error="0.0,0.0,abs"
+                             L2 norm = 0.0, Linf = inf
   --transport-parameters | -T         Specify File transport parameters
                                       e.g. "Library=stdio"
   --engine               | -E <name>  Specify ADIOS Engine