diff --git a/pyproject.toml b/pyproject.toml
index 1512bc6..d6b72d4 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "hec-dss-python"
-version = "0.1.22"
+version = "0.1.23"
 description = "Python wrapper for the HEC-DSS file database C library."
 authors = ["Hydrologic Engineering Center"]
 license = "MIT"
diff --git a/src/hecdss/hecdss.py b/src/hecdss/hecdss.py
index 16cdfe2..164ab60 100644
--- a/src/hecdss/hecdss.py
+++ b/src/hecdss/hecdss.py
@@ -301,6 +301,13 @@ def _get_paired_data(self, pathname):
         pd.ordinates = np.array(doubleOrdinates)
 
         n = numberCurves2[0].value
+        ordinateCount = len(doubleOrdinates)
+        if n > 1:
+            # ---------------------------------------------------------------------------------------- #
+            # rearrange from consecutive values for each curve to consecutive curves for each ordinate #
+            # ---------------------------------------------------------------------------------------- #
+            groups = [doubleValues[i*ordinateCount:(i+1)*ordinateCount] for i in range(n)]
+            doubleValues = list(map(list, zip(*groups)))
         pd.values = np.array(doubleValues).reshape((len(doubleOrdinates), n))
         # pd.values = [doubleValues[i:i+n] for i in range(0, len(doubleValues), n)]
         pd.labels = labels
@@ -515,7 +522,7 @@ def _get_location_info(self, pathname: str):
         )
 
         if status != 0:
-            print("Function call failed with result:", status)
+            # print("Function call failed with result:", status)
             return None
 
         location_info = LocationInfo.create(
diff --git a/src/hecdss/native.py b/src/hecdss/native.py
index 61ce2e0..1a0ceff 100644
--- a/src/hecdss/native.py
+++ b/src/hecdss/native.py
@@ -559,14 +559,22 @@ def hec_dss_pdStore(
             c_char_p,  # timeZoneName (const char*) - New argument
         ]
 
+        numberCurves = len(pd.values[0])
+        if numberCurves > 1:
+            _values = pd.values.tolist()
+            if len(_values[0]) > 1:
+                _values = [[_values[i][j] for i in range(len(_values))] for j in range(len(_values[0]))]
+            values2 = np.array(_values)
+        else:
+            values2 = pd.values
         c_pathname = c_char_p(pd.id.encode("utf-8"))
         c_Ordinates = (c_double * len(pd.ordinates))(*pd.ordinates)
         c_OrdinatesLength = len(pd.ordinates)
-        flat_list = pd.values.flatten()
+        flat_list = values2.flatten()
         c_Values = (c_double * len(flat_list))(*flat_list)
         c_ValuesLength = len(flat_list)
         c_numberOrdinates = len(pd.ordinates)
-        c_numberCurves = len(pd.values[0])
+        c_numberCurves = numberCurves
         c_unitsIndependent = c_char_p(pd.units_independent.encode("utf-8"))
         c_typeIndependent = c_char_p(pd.type_independent.encode("utf-8"))
         c_unitsDependent = c_char_p(pd.units_dependent.encode("utf-8"))
@@ -1094,8 +1102,8 @@ def hec_dss_locationRetrieve(self, fullPath: str, x: List[float], y: List[float]
             verticalDatum[0] = c_verticalDatum.value
             timeZoneName[0] = c_timeZoneName.value.decode("utf-8")
             supplemental[0] = c_supplemental.value.decode("utf-8")
-        else:
-            print("Unable to read location information, Function call failed with result:", result)
+        # else:
+        #     print("Unable to read location information, Function call failed with result:", result)
 
         return result