diff --git a/src/Solver/Mage.py b/src/Solver/Mage.py
index 56dca89a8c32eb0cfbc8ce1f52a0f0209097fbae..75a03949bd230fbbf266c7494c1bd14c0f2f4357 100644
--- a/src/Solver/Mage.py
+++ b/src/Solver/Mage.py
@@ -998,7 +998,7 @@ class Mage8(Mage):
             logger.debug(f"read_gra: nb_profile = {nb_profile}")
             logger.debug(f"read_gra: mage_version = {mage_version}")
 
-            if mage_version <= 80:
+            if mage_version < 80:
                 msg = (
                     "Read GRA files: " +
                     f"Possible incompatible mage version '{mage_version}', " +
@@ -1066,60 +1066,69 @@ class Mage8(Mage):
             def ip_to_ri(r, i): return i - reach_offset[r]
 
             ts = set()
-            ind = 0
             end = False
 
             newline()
             while not end:
                 n = read_int(1)[0]
                 timestamp = read_float64(1)[0]
+                with_bedload = read_int(1)[0]
 
-                logger.debug(f"read_gra: timestamp = {timestamp} sec")
-                ts.add(timestamp)
+                logger.debug(f"read_gra: Number of cross section: {n}")
+                logger.debug(f"read_gra: Timestamp: {timestamp}")
+                logger.debug(f"read_gra: Type of bedload: {with_bedload}")
 
                 endline()
 
-                for i in range(n):
+                npts = [1] * n
+                if with_bedload == 1:
                     newline()
-                    nsl = read_int(1)[0]
+                    npts = read_int(n)
                     endline()
+                sum_npts = sum(npts)
+                logger.debug(f"read_gra: Number of points: {sum_npts}")
 
-                    # Get current profile id
-                    reach = ip_to_r(i)
-                    ri = ip_to_ri(reach, i)
+                newline()
+                nsl = read_int(sum_npts)
+                logger.debug(f"read_gra: Number of sedimentary layers: {nsl}")
+                endline()
 
-                    if nsl > 1:
-                        logger.warning(
-                            "read_gra: " +
-                            "Multiple sediment layers for one profile " +
-                            "is not implemented yet..."
-                        )
+                newline()
+                data = read_float64(sum(nsl) * 3)
+                endline()
 
-                    for j in range(nsl):
-                        newline()
-                        nl = read_int(1)[0]
-                        endline()
+                ts.add(timestamp)
 
-                        sl = []
+                i_pts = 0
+                i_data = 0
+                # Loop on cross section
+                for i in range(n):
+                    sec_sl = []
+                    reach = ip_to_r(i)
+                    p_i = ip_to_ri(reach, i)
 
-                        for k in range(nl):
-                            newline()
-                            data = read_float64(3)
-                            endline()
+                    # Loop on cross section points
+                    for j in range(npts[i]):
+                        sl = []
 
-                            h = data[0]
-                            d50 = data[1]
-                            sigma = data[2]
+                        # Loop on sediment layers
+                        for k in range(nsl[i_pts]):
+                            h = data[i_data]
+                            d50 = data[i_data + 1]
+                            sigma = data[i_data + 2]
 
                             sl.append((h, d50, sigma))
+                            i_data += 3
 
-                        reach.set(ri, timestamp, "sl", sl)
+                        i_pts += 1
+                        sec_sl.append(sl)
 
-                ind += 1
+                    reach.set(p_i, timestamp, "sl", sec_sl)
+
+                logger.debug(f"read_gra: data size = {len(data)} ({i_data} readed)")
                 end = newline().size <= 0
 
-            logger.debug(reachs[0].profiles[0]._data)
-            results.set("timestamps", ts)
+            results.set("sediment_timestamps", ts)
             logger.info(f"read_gra: ... end with {len(ts)} timestamp read")
 
     @timer