From 60346b3ebe92fd07dd6ef9c242319934509c27de Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Tue, 21 Apr 2026 18:20:36 +0200
Subject: [PATCH 01/16] [no ci] WIP name may change add metafspm model class

---
 src/openalea/hydroroot/main.py | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index c475985..125417b 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -11,6 +11,27 @@
 from openalea.hydroroot.generator import markov, measured_root # 21-12-14: FB __init__.py in src not doing job
 from openalea.hydroroot.water_solute_transport import pressure_calculation_no_non_permeating_solutes, init_some_MTG_properties, pressure_calculation
 
+from dataclasses import dataclass, field, fields
+
+from openalea.metafspm.component import Model, declare
+from openalea.metafspm.component_factory import *
+
+@dataclass
+class hydroroot_model(Model):
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+
 def hydroroot_mtg(
     primary_length=0.15,
     delta=2.e-3,

From 8045725db920acbd9f0ede8d0ce343bf14c3d0e0 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Tue, 21 Apr 2026 18:21:23 +0200
Subject: [PATCH 02/16] [no ci] remove volume calculation in root_builder not
 returned and not used

---
 src/openalea/hydroroot/main.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index 125417b..877fa53 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -486,7 +486,7 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     # total_length is the total length of the RSA (sum of the length of all the segments)
     total_length = g.nb_vertices(scale = 1) * segment_length
     g, surface = radius.compute_surface(g)
-    g, volume = radius.compute_volume(g)
+    # g, volume = radius.compute_volume(g)
 
     if df is not None:
         v_base = next(g.component_roots_at_scale_iter(g.root, scale = g.max_scale()))

From ad9127d003d061311846d20910fd44549bf23bdd Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 09:34:03 +0200
Subject: [PATCH 03/16] [no ci] avoid specific column names in length law issue
 #110

---
 src/openalea/hydroroot/generator/markov.py | 2 +-
 src/openalea/hydroroot/law.py              | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/openalea/hydroroot/generator/markov.py b/src/openalea/hydroroot/generator/markov.py
index 249bfdb..8a18b3d 100644
--- a/src/openalea/hydroroot/generator/markov.py
+++ b/src/openalea/hydroroot/generator/markov.py
@@ -309,7 +309,7 @@ def generate_g(seed = None, length_data = None, branching_variability = 0.25,
 
     # F. Bauget 2022-08-12: added if-else to be able to use the function without length data, useful for usage demo
     if length_data:
-        length_max_secondary = length_data[0].LR_length_mm.max() * 1e-3  # in m
+        length_max_secondary = length_data[0].iloc[:,0].max() * 1e-3  # in m
 
         law_order1 = length_law(length_data[0], scale_x = primary_length / 100., scale = segment_length)
         law_order2 = length_law(length_data[1], scale_x = length_max_secondary / 100., scale = segment_length)
diff --git a/src/openalea/hydroroot/law.py b/src/openalea/hydroroot/law.py
index 3401243..094c39d 100644
--- a/src/openalea/hydroroot/law.py
+++ b/src/openalea/hydroroot/law.py
@@ -252,8 +252,8 @@ def length_law(pd, scale_x = 1 / 100., scale_y = 1., scale = 1e-4, uniform = 'ex
     * 2nd col: "relative_distance_to_tip" relative distance to tip in % so between 0 and 100.
 
     """
-    x = pd.relative_distance_to_tip.tolist()
-    y = pd.LR_length_mm.tolist()
+    x = pd.iloc[:,1].tolist()
+    y = pd.iloc[:,0].tolist()
 
     # size of the windows: in %
     size *= scale_x

From a55e16245b74d6b41650fc64befb5444d5b200d3 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 09:37:36 +0200
Subject: [PATCH 04/16] [no ci] fix a todo in length_law

---
 src/openalea/hydroroot/law.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/openalea/hydroroot/law.py b/src/openalea/hydroroot/law.py
index 094c39d..f41722b 100644
--- a/src/openalea/hydroroot/law.py
+++ b/src/openalea/hydroroot/law.py
@@ -232,7 +232,7 @@ def reference_relative_law(x, y, size=5e-2, scale_x=1., scale_y=1e-3):
     return length.fit_law(X, means, ext=2)
 
 
-def length_law(pd, scale_x = 1 / 100., scale_y = 1., scale = 1e-4, uniform = 'expo', size = 5):
+def length_law(pd, scale_x = 1 / 100., scale_y = 1.e-3, scale = 1e-4, uniform = 'expo', size = 5):
     """Build the function giving the lateral length according to its position on the parent branch
 
     :param pd: DataFrame
@@ -257,11 +257,11 @@ def length_law(pd, scale_x = 1 / 100., scale_y = 1., scale = 1e-4, uniform = 'ex
 
     # size of the windows: in %
     size *= scale_x
-    # TODO : change '1.e-3 * scale_y' to scale_y
+    
     _length_law = histo_relative_law(x, y,
                                      size = size,
                                      scale_x = scale_x,
-                                     scale_y = 1.e-3 * scale_y,
+                                     scale_y = scale_y,
                                      scale = scale,
                                      plot = False,
                                      uniform = uniform)

From 341409bee5346f531f7db2429f17614a60b1d905 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 09:43:05 +0200
Subject: [PATCH 05/16] [no ci] added a todo

---
 src/openalea/hydroroot/generator/markov.py | 1 +
 src/openalea/hydroroot/law.py              | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/openalea/hydroroot/generator/markov.py b/src/openalea/hydroroot/generator/markov.py
index 8a18b3d..7865f8d 100644
--- a/src/openalea/hydroroot/generator/markov.py
+++ b/src/openalea/hydroroot/generator/markov.py
@@ -309,6 +309,7 @@ def generate_g(seed = None, length_data = None, branching_variability = 0.25,
 
     # F. Bauget 2022-08-12: added if-else to be able to use the function without length data, useful for usage demo
     if length_data:
+        #TODO avoid unit management here mm to m '* 1e-3'
         length_max_secondary = length_data[0].iloc[:,0].max() * 1e-3  # in m
 
         law_order1 = length_law(length_data[0], scale_x = primary_length / 100., scale = segment_length)
diff --git a/src/openalea/hydroroot/law.py b/src/openalea/hydroroot/law.py
index f41722b..1c7101e 100644
--- a/src/openalea/hydroroot/law.py
+++ b/src/openalea/hydroroot/law.py
@@ -257,7 +257,7 @@ def length_law(pd, scale_x = 1 / 100., scale_y = 1.e-3, scale = 1e-4, uniform =
 
     # size of the windows: in %
     size *= scale_x
-    
+
     _length_law = histo_relative_law(x, y,
                                      size = size,
                                      scale_x = scale_x,

From 93db6dedbb415371eb98792a76fefb6fb0736053 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 10:40:37 +0200
Subject: [PATCH 06/16] [no ci] fix docstring

---
 src/openalea/hydroroot/main.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index 877fa53..3bf0047 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -437,7 +437,7 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     :param nude_length: (float) - length from tip without lateral for generated mtg (Default value = 2.0e-2)
     :param df: (DataFrame) - DataFrame with the architecture data to be reconstructed if not None (Default value = None)
     :param segment_length: (float) - vertices length in hydroroot should be in m (Default value = 1.0e-4)
-    :param length_data: (string list) - the file name with the length data laws (Default value = None)
+    :param length_data: (list of DataFrame) - DataFrame with the length data laws (Default value = None)
     :param branching_variability: (float) (Default value = 0.25)
     :param order_max: (float) (Default value = 4.0)
     :param order_decrease_factor: (float) (Default value = 0.7)

From 2675992cc3715bf49d0c109f33234478a85a1e8a Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 14:45:32 +0200
Subject: [PATCH 07/16] [no ci] moved model class to from main.py to model.py

---
 src/openalea/hydroroot/main.py  |  21 ------
 src/openalea/hydroroot/model.py | 112 ++++++++++++++++++++++++++++++++
 2 files changed, 112 insertions(+), 21 deletions(-)
 create mode 100644 src/openalea/hydroroot/model.py

diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index 3bf0047..dc68296 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -11,27 +11,6 @@
 from openalea.hydroroot.generator import markov, measured_root # 21-12-14: FB __init__.py in src not doing job
 from openalea.hydroroot.water_solute_transport import pressure_calculation_no_non_permeating_solutes, init_some_MTG_properties, pressure_calculation
 
-from dataclasses import dataclass, field, fields
-
-from openalea.metafspm.component import Model, declare
-from openalea.metafspm.component_factory import *
-
-@dataclass
-class hydroroot_model(Model):
-
-    def __init__(self, g, time_step, **scenario):
-        # 4/21/26 FB: copied from root_cynaps RootWaterModel
-        self.g = g
-        self.props = self.g.properties()
-        self.time_step = time_step
-        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
-        self.vertices = self.g.vertices(scale=self.g.max_scale())
-
-        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
-        self.apply_scenario(**scenario)
-        self.link_self_to_mtg()
-
-
 def hydroroot_mtg(
     primary_length=0.15,
     delta=2.e-3,
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
new file mode 100644
index 0000000..47049d1
--- /dev/null
+++ b/src/openalea/hydroroot/model.py
@@ -0,0 +1,112 @@
+from dataclasses import dataclass
+
+from openalea.metafspm.component import Model, declare
+
+from openalea.hydroroot.law import length_law
+
+@dataclass
+class RsaModel(Model):
+    # Input Parameters related to architecture
+    length_data: list = declare(default=None, unit="[adim,m]", unit_comment="relative distance to the tip, meter",
+                                description="list of float or list of 2 list of float, lateral length vs relat. dist. to tip",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    primary_length: float = declare(default=0.16, unit="m", unit_comment="", description="length of the primary root",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    seed: int = declare(default=None, unit="", unit_comment="", description="seed used to generate the architecture",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    branching_delay: float = declare(default=2.0e-3, unit="m", unit_comment="",
+                                     description="distance between branching laterals",
+                                     min_value="0", max_value="", value_comment="", references="", DOI="",
+                                     variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    branching_variability: float = declare(default=0.25, unit="", unit_comment="",
+                                           description="variability of the branching laterals distance",
+                                           min_value="0", max_value="1", value_comment="", references="", DOI="",
+                                           variable_type="parameter", by="RsaModel", state_variable_type="",
+                                           edit_by="user")
+
+    order_max: int = declare(default=4, unit="", unit_comment="", description="maximum order of laterals",
+                             min_value="0", max_value="1", value_comment="", references="", DOI="",
+                             variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    segment_length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    nude_length: float = declare(default=0.02, unit="m", unit_comment="",
+                                 description="part of roots without any lateral root, distance from tip",
+                                 min_value="0", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    ref_radius: float = declare(default=7.0e-5, unit="m", unit_comment="",
+                                description="reference radius or radius of the primary root",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    order_decrease_factor: float = declare(default=0.7, unit="", unit_comment="",
+                                           description="radius decrease factor applied when increasing order",
+                                           min_value="0", max_value="1", value_comment="", references="", DOI="",
+                                           variable_type="parameter", by="RsaModel", state_variable_type="",
+                                           edit_by="user")
+
+    # Plant scale variable
+    surface: float = declare(default=None, unit="m2", unit_comment="", description="total surface of the RSA",
+                             min_value="0", max_value="", value_comment="", references="", DOI="",
+                             variable_type="plant_scale_state", by="RsaModel", state_variable_type="", edit_by="user")
+
+    volume: float = declare(default=None, unit="m3", unit_comment="", description="total surface of the RSA",
+                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                            variable_type="plant_scale_state", by="RsaModel", state_variable_type="", edit_by="user")
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+    def generate_mtg(self):
+        """
+        generate a MTG according to the input parameters using self.props
+        """
+        nb_vertices = int(self.primary_length / self.segment_length)
+        # some distance expressed in term of nb of vertices
+        nb_nude_vertices = int(self.nude_length / self.segment_length)
+        nb_branching_delay = int(self.branching_delay / self.segment_length)
+
+        if isinstance(self.length_data, list):
+            length_max_secondary = self.length_data[0].iloc[:, 0]
+            law_order1 = length_law(self.length_data[0], scale_x=self.primary_length / 100., scale=self.segment_length)
+            law_order2 = length_law(self.length_data[1], scale_x=length_max_secondary / 100., scale=self.segment_length)
+            _length_law = [law_order1, law_order2]
+
+        self.g = markov.markov_binary_tree(nb_vertices=nb_vertices,
+                                           branching_variability=self.branching_variability,
+                                           branching_delay=nb_branching_delay,
+                                           length_law=length_law,
+                                           nude_tip_length=nb_nude_vertices,
+                                           order_max=self.order_max,
+                                           seed=self.seed)
+        _g, self.surface = radius.compute_surface(self.g)
+        _g, self.volume = radius.compute_volume(self.g)
+
+    def set_mtg(self, _g):
+        """
+        set self.g and compute surface, volume, etc.
+        :param _g: (MTG) existing mtg, for instance, comming from file reading
+        """
+        self.g = _g
+        _g, self.surface = radius.compute_surface(self.g)
+        _g, self.volume = radius.compute_volume(self.g)
+

From da5601607a94471a78789607fb197d88fb54e49e Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 15:22:52 +0200
Subject: [PATCH 08/16] [no ci] fix import in model.py

---
 src/openalea/hydroroot/model.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index 47049d1..b351799 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -2,6 +2,8 @@
 
 from openalea.metafspm.component import Model, declare
 
+from openalea.hydroroot.generator import markov
+from openalea.hydroroot import radius
 from openalea.hydroroot.law import length_law
 
 @dataclass

From 6ed4f73996fef4ec6befd32c41c27d6f52db6558 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 15:44:22 +0200
Subject: [PATCH 09/16] [no ci] fix issue #111

changed property name mylength to dist_to_base
---
 doc/example/example_rsml_io.py           | 10 +++++-----
 src/openalea/hydroroot/analysis.py       |  6 +++---
 src/openalea/hydroroot/display.py        |  8 ++++----
 src/openalea/hydroroot/flux.py           |  6 +++---
 src/openalea/hydroroot/main.py           |  4 ++--
 src/openalea/hydroroot/solver_wrapper.py |  4 ++--
 test/test_rsml_yaml.py                   | 10 +++++-----
 7 files changed, 24 insertions(+), 24 deletions(-)

diff --git a/doc/example/example_rsml_io.py b/doc/example/example_rsml_io.py
index 0677e80..3f2897e 100644
--- a/doc/example/example_rsml_io.py
+++ b/doc/example/example_rsml_io.py
@@ -13,11 +13,11 @@
   into continuous mtg representation
 - set the resolution according to the unit precised in the rsml, in Hydroroot lengths are in meter
 - transform the continuous mtg to discrete mtg usable in Hydroroot according to resolution
-- calculation of g properties (radius, mylength) and  flux
+- calculation of g properties (radius, dist_to_base) and  flux
 
 - export discrete mtg to rsml
 - re-import rsml to continuous mtg, transform the latter to discrete mtg
-- re-do calculation of g properties (radius, mylength) and  flux
+- re-do calculation of g properties (radius, dist_to_base) and  flux
 
 - compare both calculations should be zero
 
@@ -53,7 +53,7 @@ def set_mtg_properties(g):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -79,7 +79,7 @@ def set_mtg_properties(g):
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
         _mylength[v] = _mylength[pid] + parameter.archi['segment_length'] if pid else parameter.archi['segment_length']
-    g.properties()['mylength'] = _mylength
+    g.properties()['dist_to_base'] = _mylength
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * parameter.archi['segment_length']
@@ -139,7 +139,7 @@ def hydro_calculation(g, axfold = 1., radfold = 1., axial_data = None, k_radial
     # continuous mtg to discrete mtg
     g = import_rsml_to_discrete_mtg(g_c, segment_length = parameter.archi['segment_length'], resolution = resolution)
 
-    # calculation of g properties: radius, mylength, etc.
+    # calculation of g properties: radius, dist_to_base, etc.
     g, primary_length, _length, surface = set_mtg_properties(g)
 
     # flux calculation
diff --git a/src/openalea/hydroroot/analysis.py b/src/openalea/hydroroot/analysis.py
index c4718f0..0478e77 100644
--- a/src/openalea/hydroroot/analysis.py
+++ b/src/openalea/hydroroot/analysis.py
@@ -23,13 +23,13 @@ def nb_roots(g, l, root=1, dl=1e-4, max_order=None):
     """
     length = {}
 
-    if 'mylength' in g.property_names():
-        length = g.property('mylength')
+    if 'dist_to_base' in g.property_names():
+        length = g.property('dist_to_base')
     else:
         for v in pre_order2(g, root):
             pid = g.parent(v)
             length[v] = length[pid] + dl if pid else dl
-        g.properties()['mylength'] = length
+        g.properties()['dist_to_base'] = length
 
     order = None
     if max_order is not None:
diff --git a/src/openalea/hydroroot/display.py b/src/openalea/hydroroot/display.py
index 56394ac..e4c45fc 100644
--- a/src/openalea/hydroroot/display.py
+++ b/src/openalea/hydroroot/display.py
@@ -37,9 +37,9 @@ def root_visitor(g, v, turtle, prune=prune):
         :param prune: (float) - distance from base after witch the MTG is no longer read (Default value = None)
 
         """
-        mylength = {}
-        if prune and ('mylength' in g.properties()):
-            mylength = g.property('mylength')
+        dist_to_base = {}
+        if prune and ('dist_to_base' in g.properties()):
+            dist_to_base = g.property('dist_to_base')
         angles = [90,45]+[30]*5
         n = g.node(v)
         # radius = n.radius*1.e4
@@ -49,7 +49,7 @@ def root_visitor(g, v, turtle, prune=prune):
         length = n.length * factor
 
         if prune:
-            if mylength.get(v,0.)>prune:
+            if dist_to_base.get(v,0.)>prune:
                 return
 
         if g.edge_type(v) == '+':
diff --git a/src/openalea/hydroroot/flux.py b/src/openalea/hydroroot/flux.py
index fd0fb58..a70bac3 100644
--- a/src/openalea/hydroroot/flux.py
+++ b/src/openalea/hydroroot/flux.py
@@ -482,13 +482,13 @@ def segments_at_length(g, l, root=1, dl=1e-4):
     """
     length = {}
 
-    if 'mylength' in g.property_names():
-        length = g.property('mylength')
+    if 'dist_to_base' in g.property_names():
+        length = g.property('dist_to_base')
     else:
         for v in traversal.pre_order2(g, root):
             pid = g.parent(v)
             length[v] = length[pid] + dl if pid else dl
-        g.properties()['mylength'] = length
+        g.properties()['dist_to_base'] = length
 
     vids = []
     for v in g:
diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index dc68296..84a4b4a 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -407,7 +407,7 @@ def hydroroot_from_data(
 def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length = 2.0e-2, df = None, segment_length = 1.0e-4,
                   length_data = None, branching_variability = 0.25, order_max = 4.0, order_decrease_factor = 0.7,
                   ref_radius = 7.0e-5, Flag_radius = False):
-    """wrapper function: build a MTG with properties that are set like : radius, vertex length, position (distance to tip) and mylength (distance to base).
+    """wrapper function: build a MTG with properties that are set like : radius, vertex length, position (distance to tip) and dist_to_base (distance to base).
     The MTG is either generated or built from a DataFrame. The unit length in hydroroot should be in m.
 
     :param primary_length: (float) - primary root length for generated mtg (Default value = 0.13)
@@ -460,7 +460,7 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
         _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+    g.properties()['dist_to_base'] = _mylength
 
     # total_length is the total length of the RSA (sum of the length of all the segments)
     total_length = g.nb_vertices(scale = 1) * segment_length
diff --git a/src/openalea/hydroroot/solver_wrapper.py b/src/openalea/hydroroot/solver_wrapper.py
index 4434ef8..c5cfcef 100644
--- a/src/openalea/hydroroot/solver_wrapper.py
+++ b/src/openalea/hydroroot/solver_wrapper.py
@@ -731,8 +731,8 @@ def Jv_cnf_calculation(g, sigma, Js, Ps):
 
         # case where the primary is shorter than laterals
         max_length = primary_length
-        mylength = g_cut[0, ig].property('mylength')
-        if max(mylength.values()) > max_length: max_length = max(mylength.values())
+        dist_to_base = g_cut[0, ig].property('dist_to_base')
+        if max(dist_to_base.values()) > max_length: max_length = max(dist_to_base.values())
 
         # set conductance
         g_cut[0, ig] = set_conductances(g_cut[0, ig], axial_pr = axial_data, k0_pr = k[0], axial_lr = axial_lr,
diff --git a/test/test_rsml_yaml.py b/test/test_rsml_yaml.py
index 949f64b..183721c 100644
--- a/test/test_rsml_yaml.py
+++ b/test/test_rsml_yaml.py
@@ -18,7 +18,7 @@ def root_creation(g, segment_length, ref_radius, order_decrease_factor):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -44,7 +44,7 @@ def root_creation(g, segment_length, ref_radius, order_decrease_factor):
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
         _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+    g.properties()['dist_to_base'] = _mylength
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length
@@ -82,7 +82,7 @@ def test_rsml():
     resolution *= rsml_units_to_metre[unit]  # rsml file unit to meter
     g = import_rsml_to_discrete_mtg(g_c, segment_length = segment_length, resolution = resolution)
     
-    # calculation of g properties: radius, mylength, etc.
+    # calculation of g properties: radius, dist_to_base, etc.
     g, primary_length, _length, surface = set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor)
 
     export_mtg_to_rsml(g, "data/test_rsml_io.rsml", segment_length = segment_length)
@@ -135,7 +135,7 @@ def set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor):
     Set MTG properties and perform some gemetrical calculation
 
     The vertex radius properties is set.
-    The following properties are computed: length, position, mylength, surface, volume, total length,
+    The following properties are computed: length, position, dist_to_base, surface, volume, total length,
         primary root length
 
     :param:
@@ -161,7 +161,7 @@ def set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor):
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
         _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['mylength'] = _mylength
+    g.properties()['dist_to_base'] = _mylength
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length

From 763760537fda39ffe1a3d52b278c667bd7e36087 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 15:49:34 +0200
Subject: [PATCH 10/16] [no ci] gitignore

---
 .gitignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 91cc5cb..60a421c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,4 +77,4 @@ doc/_build
 
 # baugetfa
 tmp*
-
+test/data/test_rsml_io.rsml

From 81c8ee0b917ef8c1675e2f61c44287a3b4bf5fbd Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 16:40:49 +0200
Subject: [PATCH 11/16] [no ci] changed local var name mylength to dist_to_base

---
 doc/example/example_rsml_io.py   |  6 ++--
 src/openalea/hydroroot/main.py   |  6 ++--
 src/openalea/hydroroot/model.py  | 53 +++++++++++++++++++++++++++++---
 src/openalea/hydroroot/radius.py |  3 ++
 test/test_rsml_yaml.py           | 12 ++++----
 5 files changed, 63 insertions(+), 17 deletions(-)

diff --git a/doc/example/example_rsml_io.py b/doc/example/example_rsml_io.py
index 3f2897e..b86879b 100644
--- a/doc/example/example_rsml_io.py
+++ b/doc/example/example_rsml_io.py
@@ -75,11 +75,11 @@ def set_mtg_properties(g):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + parameter.archi['segment_length'] if pid else parameter.archi['segment_length']
-    g.properties()['dist_to_base'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + parameter.archi['segment_length'] if pid else parameter.archi['segment_length']
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * parameter.archi['segment_length']
diff --git a/src/openalea/hydroroot/main.py b/src/openalea/hydroroot/main.py
index 84a4b4a..cdba8aa 100644
--- a/src/openalea/hydroroot/main.py
+++ b/src/openalea/hydroroot/main.py
@@ -456,11 +456,11 @@ def root_builder(primary_length = 0.13, seed = None, delta = 2.0e-3, nude_length
     g = radius.compute_relative_position(g)
 
     # Calculation of the distance from base of each vertex, used for cut and flow
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['dist_to_base'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # total_length is the total length of the RSA (sum of the length of all the segments)
     total_length = g.nb_vertices(scale = 1) * segment_length
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index b351799..579d5cd 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -100,15 +100,58 @@ def generate_mtg(self):
                                            nude_tip_length=nb_nude_vertices,
                                            order_max=self.order_max,
                                            seed=self.seed)
-        _g, self.surface = radius.compute_surface(self.g)
-        _g, self.volume = radius.compute_volume(self.g)
 
-    def set_mtg(self, _g):
+        self.compute_surface_volume()
+        self.compute_radius()
+
+    def compute_length(self):
+        """
+        compute vertices length according to self.segment_length
+        :return:
+        """
+        self.g = radius.compute_length(self.g, self.segment_length)
+
+    def compute_dist_to_base(self):
+        """
+
+        :return:
+        """
+
+    def compute_dist_to_tip(self):
+        """
+        For each vertex compute the distance from the tip of the axis bearing it in absolute
+        and relative to the length of the axis bearing it
+        :return:
+        """
+        self.g = radius.compute_relative_position(self.g)
+
+    def compute_radius(self):
+        """
+        compute the radius of each vertex from self.ref_radius with fixed decrease between each order.
+        :return:
+        """
+        self.g = radius.ordered_radius(self.g, self.ref_radius, self.order_decrease_factor)
+
+    def compute_surface_volume(self):
         """
         set self.g and compute surface, volume, etc.
-        :param _g: (MTG) existing mtg, for instance, comming from file reading
         """
-        self.g = _g
+
         _g, self.surface = radius.compute_surface(self.g)
         _g, self.volume = radius.compute_volume(self.g)
 
+@dataclass
+class HydrostaticModel(Model):
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
diff --git a/src/openalea/hydroroot/radius.py b/src/openalea/hydroroot/radius.py
index fcf66f2..fe21792 100644
--- a/src/openalea/hydroroot/radius.py
+++ b/src/openalea/hydroroot/radius.py
@@ -223,6 +223,9 @@ def compute_relative_position(g):
     return:
         - g
 
+    .. note::
+        For instance, needed to compute axial and radial conductances that are properties depending
+        on the distance to tip
     """
     #print 'entering MTG node positionning computation'
     scale = g.max_scale()
diff --git a/test/test_rsml_yaml.py b/test/test_rsml_yaml.py
index 183721c..18149a2 100644
--- a/test/test_rsml_yaml.py
+++ b/test/test_rsml_yaml.py
@@ -40,11 +40,11 @@ def root_creation(g, segment_length, ref_radius, order_decrease_factor):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['dist_to_base'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length
@@ -157,11 +157,11 @@ def set_mtg_properties(g, segment_length, ref_radius, order_decrease_factor):
     # Calculation of the distance from base of each vertex, used for cut and flow
     # Remark: this calculation is done in flux.segments_at_length; analysis.nb_roots but there is a concern with the
     # parameter dl which should be equal to vertex length but which is not pass
-    _mylength = {}
+    _dist_to_base = {}
     for v in traversal.pre_order2(g, 1):
         pid = g.parent(v)
-        _mylength[v] = _mylength[pid] + segment_length if pid else segment_length
-    g.properties()['dist_to_base'] = _mylength
+        _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+    g.properties()['dist_to_base'] = _dist_to_base
 
     # _length is the total length of the RSA (sum of the length of all the segments)
     _length = g.nb_vertices(scale = 1) * segment_length

From 84f3d70d092909249595a1b9e97afcee5319c2fd Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Wed, 22 Apr 2026 17:59:12 +0200
Subject: [PATCH 12/16] [no ci] HydrostaticModel and Soil1DModel

---
 .gitignore                      |   2 +
 src/openalea/hydroroot/model.py | 151 +++++++++++++++++++++++++++++++-
 2 files changed, 149 insertions(+), 4 deletions(-)

diff --git a/.gitignore b/.gitignore
index 60a421c..380068f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -78,3 +78,5 @@ doc/_build
 # baugetfa
 tmp*
 test/data/test_rsml_io.rsml
+out.csv
+
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index 579d5cd..443974c 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -1,10 +1,13 @@
 from dataclasses import dataclass
 
+from openalea.mtg import traversal
+
 from openalea.metafspm.component import Model, declare
 
+from openalea.hydroroot import radius, conductance, main, flux
 from openalea.hydroroot.generator import markov
-from openalea.hydroroot import radius
 from openalea.hydroroot.law import length_law
+from openalea.hydroroot.length import fit_law
 
 @dataclass
 class RsaModel(Model):
@@ -57,6 +60,19 @@ class RsaModel(Model):
                                            variable_type="parameter", by="RsaModel", state_variable_type="",
                                            edit_by="user")
 
+    # hydraulic parameters
+    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
+    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]",
+                                            unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="parameter", by="RsaModel", state_variable_type="",
+                                            edit_by="user")
+
     # Plant scale variable
     surface: float = declare(default=None, unit="m2", unit_comment="", description="total surface of the RSA",
                              min_value="0", max_value="", value_comment="", references="", DOI="",
@@ -101,6 +117,9 @@ def generate_mtg(self):
                                            order_max=self.order_max,
                                            seed=self.seed)
 
+        self.compute_radius()
+        self.compute_length()
+        self.compute_dist_to_tip()
         self.compute_surface_volume()
         self.compute_radius()
 
@@ -113,9 +132,16 @@ def compute_length(self):
 
     def compute_dist_to_base(self):
         """
-
+        Calculation of the distance from base of each vertex, used for cut and flow
         :return:
         """
+        g = self.g
+        _dist_to_base = {}
+        for v in traversal.pre_order2(g, 1):
+            pid = g.parent(v)
+            segment_length = g.length(pid)
+            _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
+        g.properties()['dist_to_base'] = _dist_to_base
 
     def compute_dist_to_tip(self):
         """
@@ -136,12 +162,123 @@ def compute_surface_volume(self):
         """
         set self.g and compute surface, volume, etc.
         """
+        self.g, self.surface = radius.compute_surface(self.g)
+        self.g, self.volume = radius.compute_volume(self.g)
 
-        _g, self.surface = radius.compute_surface(self.g)
-        _g, self.volume = radius.compute_volume(self.g)
+    def compute_axial_conductance(self):
+        """
+        compute axial conductance of each vertex from self.ref_axial_conductivity
+        :return:
+        """
+        if self.axial_conductivity_data:
+            # Compute K using axial conductance data
+            xa, ya = self.axial_conductivity_data
+            axial_conductivity_law = fit_law(xa, ya)
+
+            self.g = conductance.fit_property_from_spline(self.g, axial_conductivity_law, 'position', 'K_exp')
+            self.g = conductance.compute_K(self.g)
+
+    def compute_radial_conductance(self):
+        """
+        compute radial conductance of each vertex from self.ref_radial_conductivity
+        :return:
+        """
+        if self.radial_conductivity_data:
+            xr, yr = self.radial_conductivity_data
+            radial_conductivity_law = fit_law(xr, yr)
+
+            self.g = conductance.fit_property_from_spline(self.g, radial_conductivity_law, 'position', 'k0')
+            self.g = conductance.compute_k(self.g, k0='k0')
 
 @dataclass
 class HydrostaticModel(Model):
+    # Hydraulic input from RsaModel
+    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
+
+    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="input", by="RsaModel", state_variable_type="",
+                                            edit_by="user")
+
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Homogeneous hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # HydrostaticModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                                description="Hydrostatic potential at the root base",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                                description="water flux at the root base, input if invert_model=False ",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    invert_model: bool = declare(default=True, unit="", unit_comment="",
+                                description="when false, distribute output flux within the root ; "
+                                            "when true, compute the output flux for the given root and conditions.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                description="Use specific model to compute conductance at tips with cut & flow.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+
+
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+    def flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+
+        Keqs = self.g.property('Keq')
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        if _psi_e is None:
+            peq = self.g.property('Peq')
+            _psi_e = peq[v_base]
+
+        self.Jv = Keqs[v_base] * (_psi_e - self.psi_base)
+
+
+@dataclass
+class Soil1DModel(Model):
+    # Soil1DModel parameters
+    soil_data: tuple = declare(default=None, unit="[m, MPa]", unit_comment="",
+                                description="tuple of 2 lists, (z,psi_e) z=depth, psi_e=water potential",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # state variable
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Soil hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="state_variable", by="HydrostaticModel", state_variable_type="", edit_by="user")
 
     def __init__(self, g, time_step, **scenario):
         # 4/21/26 FB: copied from root_cynaps RootWaterModel
@@ -155,3 +292,9 @@ def __init__(self, g, time_step, **scenario):
         self.apply_scenario(**scenario)
         self.link_self_to_mtg()
 
+    def compute_doil_psi_e(self):
+        """
+        add a soil as heterogeneous water potential along z
+        :return:
+        """
+        self.g = main.soil_1D(self.g, self.soil_data)

From d38080cf3b4f3030b30f2062a7a188e91f47af72 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Fri, 24 Apr 2026 18:46:44 +0200
Subject: [PATCH 13/16] [no ci] change models to one HydroRootModel with RSA
 hydrostatic and solute'

---
 .../example_metafspm_archi_from_file.py       | 160 +++++++++
 src/openalea/hydroroot/model.py               | 303 ++++++++++++++++--
 2 files changed, 444 insertions(+), 19 deletions(-)
 create mode 100644 doc/example/metafspm/example_metafspm_archi_from_file.py

diff --git a/doc/example/metafspm/example_metafspm_archi_from_file.py b/doc/example/metafspm/example_metafspm_archi_from_file.py
new file mode 100644
index 0000000..dbbc743
--- /dev/null
+++ b/doc/example/metafspm/example_metafspm_archi_from_file.py
@@ -0,0 +1,160 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# # Measured architecture
+
+# Read architecture from text file and run a simulation of the sap flux from an Arabidopsis de-topped root plunged in a hydroponic solution at a hydrostatic pressure of 0.4 Mpa when its  base is at the atmospheric pressure.
+
+# ## The text format architecture file
+# 
+# |distance_from_base_(mm)  |lateral_root_length_(mm)  |order|
+# |:-:  | :-:  | :-: |
+# |0.89                     |90.81             	       |1|
+# |3.02                     |63.98             	       |1|
+# |102.94                   |0.0             	         |1|
+# |2.14                     |23.72             	       |1-1|
+# |90.81                    |0.0             	         |1-1|
+# |2.48                     |5.15             	       |1-2|
+# |63.98                    |0.0             	         |1-2|
+# 
+# This is a tab separated text file with 3 columns:
+# 
+#     1. the distance from base of the branching laterals in mm
+# 
+#     2. the lateral root length in mm
+# 
+#     3. a string of one or more number indicating the parent root
+# 
+# In the example above, the root has two lateral of 1st order and on each of them one lateral of 2d order. The order of 1 indicates that the laterals are on the primary root. The last line with order 1 with 0.0 in the second column indicates the primary root tip. The line with order 1-1 indicates that this a second order lateral on the first lateral positioned at 2.14 mm from the branching on the primary root. And so on.
+# 
+# A 4th column with the averaged diameter of the root may be given, that may be used to build the MTG representing the architecture.
+
+# ## Running the calculation
+# If the package HydroRoot is not installed, the following examples can be run by cloning the sources from git and then sourcing the src directory in Ipython console for instance like this:
+
+# In[1]:
+
+
+get_ipython().run_line_magic('matplotlib', 'inline')
+
+
+# In[2]:
+
+
+import matplotlib
+from openalea.widgets.plantgl import PlantGL # notebook viewer 3D
+from openalea.plantgl.algo.view import view # 2D view
+from openalea.hydroroot.display import mtg_scene
+from openalea.hydroroot.read_file import read_archi_data
+from openalea.hydroroot.main import hydroroot_flow, root_builder
+from openalea.hydroroot.generator import measured_root
+from openalea.hydroroot import model
+from openalea.mtg import MTG
+
+
+# Read the architecture file as DataFrame
+
+# In[3]:
+
+
+df = read_archi_data('../data/plant-01.txt')
+
+
+# Building the MTG from the file, and  return the primary root length, the total length and the surface. The seed refer to the seed of the root generator when the MTG is not built from a file but is generated.
+
+# In[4]:
+
+segment_length = 1.0e-4
+time_step = 1.0e-4
+
+g = measured_root.mtg_from_aqua_data(df, segment_length)
+hydromodel = model.HydroRootModel(g, time_step)
+hydromodel.compute_metrics()
+
+# Some conductance data versus distance to tip
+
+# In[5]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[30.0,30.0]))
+hydromodel.compute_axial_conductance(data = ([0, 0.2],[3.0e-7,4.0e-4]))
+
+
+# Flux and equivalent conductance calculation, for a root in an external hydroponic medium at 0.4 MPa, its base at 0.1 MPa, and with the conductances set above.
+
+# In[6]:
+
+# g = hydromodel.g
+# g, keq, jv = hydroroot_flow(g, psi_e = 0.4, psi_base = 0.1, axial_conductivity_data = K_axial_data, radial_conductivity_data = k_radial_data)
+hydromodel.psi_e = 0.4
+hydromodel.psi_base = 0.1
+hydromodel.hydrostatic_solver_flux()
+
+# In[7]:
+
+
+print('equivalent root conductance (microL/s/MPa): ',hydromodel.keq, 'sap flux (microL/s): ', hydromodel.Jv)
+
+
+# ## Plots
+
+# Display the local water uptake heatmap in 3D
+
+# In[8]:
+
+
+s = mtg_scene(g, prop_cmap = 'j') # create a scene from the mtg with the property j is the radial flux in ul/s
+PlantGL(s) # display the root into the plantgl oawidget
+
+
+# You may change the property to display to the hydrostatic pressure inside the xylem vessels for instance
+# 
+# to reduce notebook size we use here a 2D view but you can use the openalea.widgets `PlantGL(s)` to display an interactive 3D view
+
+# In[9]:
+
+
+s = mtg_scene(g, prop_cmap='psi_in')
+view(s) #PlantGL(s)
+
+
+# You may change the radial conductivity and see the impact on the water uptake
+
+# In[10]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[300.0,300.0]))
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
+s = mtg_scene(g, prop_cmap='j')
+view(s) # to reduce notebook size we use here a 2D view but use PlantGL(s) to 3D
+
+
+# Or the axial conductance
+
+# In[11]:
+
+
+hydromodel.compute_radial_conductance(data = ([0, 0.2],[30.0,30.0]))
+hydromodel.compute_axial_conductance(data = ([0, 0.2],[3.0e-7,4.0e-4]))
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
+s = mtg_scene(g, prop_cmap='j')
+view(s) # to reduce notebook size we use here a 2D view but use PlantGL(s) to 3D
+
+
+# ## PlantGl viewer
+# 
+# This is only working on **local machine**.
+# 
+# HydroRoot has also a function using the standalone Plantgl viewer allowing to specify the mtg property.
+# 
+# The following will open a standalone window with a 3D representation of the root. Uncomment the lines to run the cell locally.
+
+# In[12]:
+
+
+# from openalea.hydroroot.display import plot
+# %gui qt
+# plot(g, prop_cmap='psi_in')
+
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index 443974c..1020dc2 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -1,16 +1,23 @@
+import math
+
 from dataclasses import dataclass
 
 from openalea.mtg import traversal
-
 from openalea.metafspm.component import Model, declare
+from openalea.metafspm.component_factory import *
 
 from openalea.hydroroot import radius, conductance, main, flux
 from openalea.hydroroot.generator import markov
 from openalea.hydroroot.law import length_law
 from openalea.hydroroot.length import fit_law
+from openalea.hydroroot.water_solute_transport import (pressure_calculation_no_non_permeating_solutes,
+                                                       init_some_MTG_properties, pressure_calculation)
+from openalea.hydroroot.analysis import intercept
+
+EPS = 1.0e-9 # epsilon value used for convergence of Newton-Raphson loop
 
 @dataclass
-class RsaModel(Model):
+class HydroRootModel(Model):
     # Input Parameters related to architecture
     length_data: list = declare(default=None, unit="[adim,m]", unit_comment="relative distance to the tip, meter",
                                 description="list of float or list of 2 list of float, lateral length vs relat. dist. to tip",
@@ -40,7 +47,7 @@ class RsaModel(Model):
                              min_value="0", max_value="1", value_comment="", references="", DOI="",
                              variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
 
-    segment_length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
+    length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
                                     min_value="0", max_value="", value_comment="", references="", DOI="",
                                     variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
 
@@ -73,7 +80,75 @@ class RsaModel(Model):
                                             variable_type="parameter", by="RsaModel", state_variable_type="",
                                             edit_by="user")
 
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                           description="Homogeneous hydrostatic potential at the vertex boundary",
+                           min_value="0", max_value="", value_comment="", references="", DOI="",
+                           variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # HydrostaticModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                              description="Hydrostatic potential at the root base",
+                              min_value="0", max_value="", value_comment="", references="", DOI="",
+                              variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                        description="water flux at the root base, input if invert_model=False ",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
+    invert_model: bool = declare(default=True, unit="", unit_comment="",
+                                 description="when false, distribute output flux within the root ; "
+                                             "when true, compute the output flux for the given root and conditions.",
+                                 min_value="", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="HydrostaticModel", state_variable_type="",
+                                 edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                 description="Use specific model to compute conductance at tips with cut & flow.",
+                                 min_value="", max_value="", value_comment="", references="", DOI="",
+                                 variable_type="parameter", by="HydrostaticModel", state_variable_type="",
+                                 edit_by="user")
+
+    # Solute parameters
+    J_s: float = declare(default=1.0e-7, unit="mol.m-2.s-1", unit_comment="", description="Active pumping rate",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    P_s: float = declare(default=1.0e-9, unit="m.s-1", unit_comment="", description="permeability coefficient",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Cse: float = declare(default=13.96, unit="mol.m-3", unit_comment="",
+                         description="concentration of permeating solutes",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Ce: float = declare(default=0.0, unit="mol.m-3", unit_comment="",
+                        description="concentration of non-permeating solutes",
+                        min_value="0", max_value="", value_comment="", references="", DOI="",
+                        variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    sigma: float = declare(default=1.0, unit="", unit_comment="", description="reflexion coefficient",
+                           min_value="0", max_value="", value_comment="", references="", DOI="",
+                           variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    temp: float = declare(default=298.0, unit="K", unit_comment="", description="sap temperature in degree K",
+                          min_value="", max_value="", value_comment="", references="", DOI="",
+                          variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    c_base: float = declare(default=None, unit="mol.m-3", unit_comment="",
+                            description="solute concentration at the root's base",
+                            min_value="", max_value="", value_comment="", references="", DOI="",
+                            variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
     # Plant scale variable
+    keq: float = declare(default=None, unit="10-9 m4.MPa-1.s-1", unit_comment="",
+                         description="equivalent conductance of the RSA",
+                         min_value="0", max_value="", value_comment="", references="", DOI="",
+                         variable_type="plant_scale_state", by="HydrostaticModel", state_variable_type="",
+                         edit_by="user")
+
     surface: float = declare(default=None, unit="m2", unit_comment="", description="total surface of the RSA",
                              min_value="0", max_value="", value_comment="", references="", DOI="",
                              variable_type="plant_scale_state", by="RsaModel", state_variable_type="", edit_by="user")
@@ -87,6 +162,7 @@ def __init__(self, g, time_step, **scenario):
         self.g = g
         self.props = self.g.properties()
         self.time_step = time_step
+        # self.choregrapher(module_family=self.__class__.__name__)
         self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
         self.vertices = self.g.vertices(scale=self.g.max_scale())
 
@@ -98,15 +174,15 @@ def generate_mtg(self):
         """
         generate a MTG according to the input parameters using self.props
         """
-        nb_vertices = int(self.primary_length / self.segment_length)
+        nb_vertices = int(self.primary_length / self.length)
         # some distance expressed in term of nb of vertices
-        nb_nude_vertices = int(self.nude_length / self.segment_length)
-        nb_branching_delay = int(self.branching_delay / self.segment_length)
+        nb_nude_vertices = int(self.nude_length / self.length)
+        nb_branching_delay = int(self.branching_delay / self.length)
 
         if isinstance(self.length_data, list):
             length_max_secondary = self.length_data[0].iloc[:, 0]
-            law_order1 = length_law(self.length_data[0], scale_x=self.primary_length / 100., scale=self.segment_length)
-            law_order2 = length_law(self.length_data[1], scale_x=length_max_secondary / 100., scale=self.segment_length)
+            law_order1 = length_law(self.length_data[0], scale_x=self.primary_length / 100., scale=self.length)
+            law_order2 = length_law(self.length_data[1], scale_x=length_max_secondary / 100., scale=self.length)
             _length_law = [law_order1, law_order2]
 
         self.g = markov.markov_binary_tree(nb_vertices=nb_vertices,
@@ -123,13 +199,31 @@ def generate_mtg(self):
         self.compute_surface_volume()
         self.compute_radius()
 
+    @actual
+    def compute_metrics(self):
+        if not self.g.property('radius'):
+            # TODO ref_radius etc.
+            self.compute_radius()
+        if not self.g.property('length'):
+            self.compute_length()
+
+        self.compute_dist_to_base()
+        self.compute_dist_to_tip()
+        self.compute_surface_volume()
+        self.g.property('surface')[1] = self.surface
+        self.g.property('volume')[1] = self.volume
+
+    @actual
+    @state
     def compute_length(self):
         """
-        compute vertices length according to self.segment_length
+        compute vertices length according to self.length
         :return:
         """
-        self.g = radius.compute_length(self.g, self.segment_length)
+        self.g = radius.compute_length(self.g, self.length)
 
+    @actual
+    @state
     def compute_dist_to_base(self):
         """
         Calculation of the distance from base of each vertex, used for cut and flow
@@ -137,12 +231,14 @@ def compute_dist_to_base(self):
         """
         g = self.g
         _dist_to_base = {}
+        segment_length = self.length
         for v in traversal.pre_order2(g, 1):
             pid = g.parent(v)
-            segment_length = g.length(pid)
             _dist_to_base[v] = _dist_to_base[pid] + segment_length if pid else segment_length
         g.properties()['dist_to_base'] = _dist_to_base
 
+    @actual
+    @state
     def compute_dist_to_tip(self):
         """
         For each vertex compute the distance from the tip of the axis bearing it in absolute
@@ -151,6 +247,8 @@ def compute_dist_to_tip(self):
         """
         self.g = radius.compute_relative_position(self.g)
 
+    @actual
+    @state
     def compute_radius(self):
         """
         compute the radius of each vertex from self.ref_radius with fixed decrease between each order.
@@ -158,6 +256,8 @@ def compute_radius(self):
         """
         self.g = radius.ordered_radius(self.g, self.ref_radius, self.order_decrease_factor)
 
+    @actual
+    @state
     def compute_surface_volume(self):
         """
         set self.g and compute surface, volume, etc.
@@ -165,11 +265,17 @@ def compute_surface_volume(self):
         self.g, self.surface = radius.compute_surface(self.g)
         self.g, self.volume = radius.compute_volume(self.g)
 
-    def compute_axial_conductance(self):
+    @actual
+    @axial
+    @state
+    def compute_axial_conductance(self, data = None):
         """
         compute axial conductance of each vertex from self.ref_axial_conductivity
         :return:
         """
+        if data:
+            self.axial_conductivity_data = data
+
         if self.axial_conductivity_data:
             # Compute K using axial conductance data
             xa, ya = self.axial_conductivity_data
@@ -178,11 +284,16 @@ def compute_axial_conductance(self):
             self.g = conductance.fit_property_from_spline(self.g, axial_conductivity_law, 'position', 'K_exp')
             self.g = conductance.compute_K(self.g)
 
-    def compute_radial_conductance(self):
+    @actual
+    @state
+    def compute_radial_conductance(self, data = None):
         """
         compute radial conductance of each vertex from self.ref_radial_conductivity
         :return:
         """
+        if data:
+            self.radial_conductivity_data = data
+
         if self.radial_conductivity_data:
             xr, yr = self.radial_conductivity_data
             radial_conductivity_law = fit_law(xr, yr)
@@ -190,9 +301,69 @@ def compute_radial_conductance(self):
             self.g = conductance.fit_property_from_spline(self.g, radial_conductivity_law, 'position', 'k0')
             self.g = conductance.compute_k(self.g, k0='k0')
 
+    def compute_intercepts(self, dists, dl=1e-4, max_order=None):
+        return intercept(self.g, dists, dl=dl, max_order=max_order)
+
+
+    @actual
+    @state
+    def hydrostatic_solver_flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.keq = self.g.property('Keq')[v_base]
+        self.Jv = self.g.property('J_out')[v_base]
+
+    @actual
+    @state
+    def solute_solver_flux(self):
+        """
+        Compute flux according to psi_e and psi_base
+        If self.psi_e is None use self.g.property('psi_e')
+        :return:
+        """
+
+        # init the MTG with no solute
+        _psi_e = self.psi_e
+        f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
+        f.run()
+        self.g = f.g
+
+        self.g = init_some_MTG_properties(self.g, tau=self.J_s, Cini=self.Cse*1.0e-9, Cpeg_ini=self.Ce*1.0e-9, t=1, Ps=self.P_s)
+
+        if self.Ce <= 0.0:
+            calculation = pressure_calculation_no_non_permeating_solutes
+        else:
+            calculation = pressure_calculation
+
+        # Newton-Raphson loop
+        nb_v = self.g.nb_vertices()
+        Fdx = 1.0
+        Fdx_old = 1.
+        while Fdx > EPS:
+            g, dx, data, row, col = calculation(g, Temp=self.temp, sigma=self.sigma, Ce=self.Ce*1.0e-9, Cse=self.Cse*1.0e-9,
+                                                Pe=self.psi_e, Pbase=self.psi_base, C_base=self.c_base*1.0e-9)
+            Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
+            if abs(Fdx - Fdx_old) < EPS: break
+            Fdx_old = Fdx
+
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.Jv = self.g.property('J_out')[v_base]
+
 @dataclass
 class HydrostaticModel(Model):
     # Hydraulic input from RsaModel
+    length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
+                                    min_value="0", max_value="", value_comment="", references="", DOI="",
+                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+
     axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
@@ -232,7 +403,87 @@ class HydrostaticModel(Model):
                                 min_value="", max_value="", value_comment="", references="", DOI="",
                                 variable_type="parameter", by="HydrostaticModel", state_variable_type="", edit_by="user")
 
+    keq: float = declare(default=None, unit="10-9 m4.MPa-1.s-1", unit_comment="", description="equivalent conductance of the RSA",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="plant_scale_state", by="HydrostaticModel", state_variable_type="", edit_by="user")
+
 
+    def __init__(self, g, time_step, **scenario):
+        # 4/21/26 FB: copied from root_cynaps RootWaterModel
+        self.g = g
+        self.props = self.g.properties()
+        self.time_step = time_step
+        self.choregrapher.add_time_and_data(instance=self, sub_time_step=self.time_step, data=self.props)
+        self.vertices = self.g.vertices(scale=self.g.max_scale())
+
+        # Before any other operation, we apply the provided scenario by changing default parameters and initialization
+        self.apply_scenario(**scenario)
+        self.link_self_to_mtg()
+
+
+@dataclass
+class SoluteModel(Model):
+    # Hydraulic input from RsaModel
+    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+                                description="(2 list of Float) axial conductivity versus dist. to tip",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
+
+    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+                                            description="(2 list of Float) radial conductivity versus dist. to tip",
+                                            min_value="0", max_value="", value_comment="", references="", DOI="",
+                                            variable_type="input", by="RsaModel", state_variable_type="",
+                                            edit_by="user")
+
+    # Soil1DModel input
+    psi_e: float = declare(default=0.401325, unit="MPa", unit_comment="",
+                                description="Homogeneous hydrostatic potential at the vertex boundary",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="input", by="Soil1DModel", state_variable_type="", edit_by="user")
+
+    # SoluteModel parameter
+    psi_base: float = declare(default=0.101325, unit="MPa", unit_comment="",
+                                description="Hydrostatic potential at the root base",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Jv: float = declare(default=None, unit="microL.s-1", unit_comment="",
+                                description="water flux at the root base, input if invert_model=False ",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    cut_and_flow: bool = declare(default=False, unit="", unit_comment="",
+                                description="Use specific model to compute conductance at tips with cut & flow.",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    J_s: float = declare(default=1.0e-7, unit="mol.m-2.s-1", unit_comment="", description="Active pumping rate",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    P_s: float = declare(default=1.0e-9, unit="m.s-1", unit_comment="", description="permeability coefficient",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Cse: float = declare(default=13.96, unit="mol.m-3", unit_comment="", description="concentration of permeating solutes",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    Ce: float = declare(default=0.0, unit="mol.m-3", unit_comment="", description="concentration of non-permeating solutes",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    sigma: float = declare(default=1.0, unit="", unit_comment="", description="reflexion coefficient",
+                                min_value="0", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    temp: float = declare(default=298.0, unit="K", unit_comment="", description="sap temperature in degree K",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
+
+    c_base: float = declare(default=None, unit="mol.m-3", unit_comment="", description="solute concentration at the root's base",
+                                min_value="", max_value="", value_comment="", references="", DOI="",
+                                variable_type="parameter", by="SoluteModel", state_variable_type="", edit_by="user")
 
     def __init__(self, g, time_step, **scenario):
         # 4/21/26 FB: copied from root_cynaps RootWaterModel
@@ -252,19 +503,33 @@ def flux(self):
         If self.psi_e is None use self.g.property('psi_e')
         :return:
         """
+
+        # init the MTG with no solute
         _psi_e = self.psi_e
         f = flux.Flux(self.g, self.Jv, _psi_e, self.psi_base, self.invert_model, cut_and_flow=self.cut_and_flow)
         f.run()
         self.g = f.g
 
-        Keqs = self.g.property('Keq')
-        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
-        if _psi_e is None:
-            peq = self.g.property('Peq')
-            _psi_e = peq[v_base]
+        self.g = init_some_MTG_properties(self.g, tau=self.J_s, Cini=self.Cse*1.0e-9, Cpeg_ini=self.Ce*1.0e-9, t=1, Ps=self.P_s)
+
+        if self.Ce <= 0.0:
+            calculation = pressure_calculation_no_non_permeating_solutes
+        else:
+            calculation = pressure_calculation
 
-        self.Jv = Keqs[v_base] * (_psi_e - self.psi_base)
+        # Newton-Raphson loop
+        nb_v = self.g.nb_vertices()
+        Fdx = 1.0
+        Fdx_old = 1.
+        while Fdx > EPS:
+            g, dx, data, row, col = calculation(g, Temp=self.temp, sigma=self.sigma, Ce=self.Ce*1.0e-9, Cse=self.Cse*1.0e-9,
+                                                Pe=self.psi_e, Pbase=self.psi_base, C_base=self.c_base*1.0e-9)
+            Fdx = math.sqrt(sum(dx ** 2.0)) / nb_v
+            if abs(Fdx - Fdx_old) < EPS: break
+            Fdx_old = Fdx
 
+        v_base = next(self.g.component_roots_at_scale_iter(self.g.root, scale=1))
+        self.Jv = self.g.property('J_out')[v_base]
 
 @dataclass
 class Soil1DModel(Model):

From aa2055e207a9473be03cf3a18b905a1122f6772c Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Fri, 24 Apr 2026 19:40:29 +0200
Subject: [PATCH 14/16] [no ci] hydroroot param to metafspm scenario

---
 .../example_metafspm_archi_from_file.py       |  10 +-
 doc/example/metafspm/parameters_plant_01.yml  | 118 ++++++++++++++++++
 src/openalea/hydroroot/init_parameter.py      |  11 ++
 src/openalea/hydroroot/model.py               |  24 ++--
 4 files changed, 150 insertions(+), 13 deletions(-)
 create mode 100644 doc/example/metafspm/parameters_plant_01.yml

diff --git a/doc/example/metafspm/example_metafspm_archi_from_file.py b/doc/example/metafspm/example_metafspm_archi_from_file.py
index dbbc743..d90e2d9 100644
--- a/doc/example/metafspm/example_metafspm_archi_from_file.py
+++ b/doc/example/metafspm/example_metafspm_archi_from_file.py
@@ -157,4 +157,12 @@
 # from openalea.hydroroot.display import plot
 # %gui qt
 # plot(g, prop_cmap='psi_in')
-
+from openalea.hydroroot.init_parameter import Parameters
+parameter = Parameters()
+parameter.read_file('parameters_plant_01.yml')
+scenario = {}
+for key, value in parameter.__dict__.items():
+    scenario.update(value)
+hydromodel = model.HydroRootModel(g, time_step, **scenario)
+hydromodel.hydrostatic_solver_flux()
+print('sap flux (microL/s): ', hydromodel.Jv)
diff --git a/doc/example/metafspm/parameters_plant_01.yml b/doc/example/metafspm/parameters_plant_01.yml
new file mode 100644
index 0000000..815f732
--- /dev/null
+++ b/doc/example/metafspm/parameters_plant_01.yml
@@ -0,0 +1,118 @@
+#%YAML
+
+#All parameters should be in SI units
+
+#Few parameters may be set to list of float or integer allowing to run successive simulation
+# there are two syntaxes:
+#       [x1, ..., xn] or range(start, end, step)
+#       eg. range(0.02, 0.09, 0.02) or [0.02, 0.04, 0.06, 0.08] will give the same results
+#       the parameter will take successively the values 0.02, 0.04, 0.06 and 0.08
+
+archi:
+    #if read_architecture is true then architecture will be constructed from the file(s) given by input_dir and input_file
+    #otherwise the architecture will be generated according to the parameters
+    read_architecture: True
+
+    #Input architecture from scanned image (distance_from_base_(mm) \t lateral_root_length_(mm) \t order)
+    #folder name
+    input_dir: ../data/
+
+    #File name:
+    #may be a list of names, eg. [file1, file2, file3] wildcar may be used
+    input_file: [plant-01.txt]
+
+    seed:
+
+    #file names with length laws relative path
+    #file format: "LR_length_mm" ; "relative_distance_to_tip"
+    #laws used to generate lateral roots of the 1st order (1_order_law), and lateral roots of order above 1 (2_order_law)
+    length_file:
+        - ../data/length*order1*.csv
+        - ../data/length*order2*.csv
+
+    #length of the primary root
+    #float or list of float
+    #unit: m
+    primary_length: 0.32
+
+    #branching delay
+    #float or list of float
+    #unit: m
+    branching_delay: 1.0e-4
+
+    #branching variability
+    #float between [0 ; 1], 0.25 means 25%
+    branching_variability: 0.25
+
+    #maximum roots order
+    order_max: 4
+
+    #vertices length
+    #unit: m
+    segment_length: 1.0e-4
+
+    #part of roots without any lateral root, distance from tip
+    #float or list of float
+    #unit: m
+    nude_length: 0.03
+
+    #reference radius of the primary root
+    #float
+    #unit: m
+    ref_radius: 7.0e-5
+
+    #radius decrease factor applied when increasing order
+    #float
+    #radius lateral order n: r = order_decrease_factor^n * ref_radius
+    order_decrease_factor: 0.7
+
+hydro:
+    #radial conductivity
+    #float
+    #unit: microL/(s.MPa.m**2)
+    k0:
+        - [0.0, 0.32]
+        - [10.0, 10.0]
+
+    #axial_conductance_data
+    #   - [x1, ......, xn]
+    #   - [K1, ....., Kn]
+    #list of float
+    #unit: microL.m/(s.Mpa)
+    axial_conductance_data:
+        - [0, 0.32]
+        - [3.0e-9,2.0e-4]
+
+experimental:
+    #water flux at the root base
+    #float
+    #unit: microL/s
+    Jv:  4.52e-3
+
+    #hydric potential outside the roots (pressure chamber)
+    #float
+    #unit: MPa
+    psi_e: 0.501325
+
+    #hydric potential at the root base (e.g. atmospheric pressure for decapitated plant)
+    #float
+    #unit: MPa
+    psi_base: .101325
+
+output:
+    #distance from the base for intercepts calculation
+    # float or list of float
+    #unit: m
+    intercepts: []
+
+    #factor to explore a k0 range
+    # float or list of float
+    radfold: 1.
+
+    #like radfold but apply to axial_conductance_data
+    axfold: 1.
+
+    #number of run with the same set of parameters i.e. number of different seeds
+    #integer
+    #enable only if read_architecture is false
+    run_nb: 1
diff --git a/src/openalea/hydroroot/init_parameter.py b/src/openalea/hydroroot/init_parameter.py
index bda715b..57da49d 100644
--- a/src/openalea/hydroroot/init_parameter.py
+++ b/src/openalea/hydroroot/init_parameter.py
@@ -200,3 +200,14 @@ def parameters_to_list(self, parameter):
                 parameter = [parameter]
 
         return parameter
+
+    def metafspm_scenario(self):
+        """
+        In metafspm scenario is just a dict
+        This function simply aggregates all variables in a single dictionary
+        :return: scenario (dict)
+        """
+        scenario = {}
+        for key, value in self.__dict__.items():
+            scenario.update(value)
+        return scenario
diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index 1020dc2..a8c881a 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -68,12 +68,12 @@ class HydroRootModel(Model):
                                            edit_by="user")
 
     # hydraulic parameters
-    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
                                 variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
 
-    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]",
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]",
                                             unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
@@ -274,11 +274,11 @@ def compute_axial_conductance(self, data = None):
         :return:
         """
         if data:
-            self.axial_conductivity_data = data
+            self.axial_conductance_data = data
 
-        if self.axial_conductivity_data:
+        if self.axial_conductance_data:
             # Compute K using axial conductance data
-            xa, ya = self.axial_conductivity_data
+            xa, ya = self.axial_conductance_data
             axial_conductivity_law = fit_law(xa, ya)
 
             self.g = conductance.fit_property_from_spline(self.g, axial_conductivity_law, 'position', 'K_exp')
@@ -292,10 +292,10 @@ def compute_radial_conductance(self, data = None):
         :return:
         """
         if data:
-            self.radial_conductivity_data = data
+            self.k0 = data
 
-        if self.radial_conductivity_data:
-            xr, yr = self.radial_conductivity_data
+        if self.k0:
+            xr, yr = self.k0
             radial_conductivity_law = fit_law(xr, yr)
 
             self.g = conductance.fit_property_from_spline(self.g, radial_conductivity_law, 'position', 'k0')
@@ -364,12 +364,12 @@ class HydrostaticModel(Model):
                                     min_value="0", max_value="", value_comment="", references="", DOI="",
                                     variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
 
-    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
                                 variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
 
-    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
                                             variable_type="input", by="RsaModel", state_variable_type="",
@@ -424,12 +424,12 @@ def __init__(self, g, time_step, **scenario):
 @dataclass
 class SoluteModel(Model):
     # Hydraulic input from RsaModel
-    axial_conductivity_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
+    axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
                                 variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
 
-    radial_conductivity_data: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
+    k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
                                             variable_type="input", by="RsaModel", state_variable_type="",

From 3dc4a52fe77806fee4c364d570e6bf0a530e64a9 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Mon, 27 Apr 2026 09:37:43 +0200
Subject: [PATCH 15/16] [no ci] fix by field

---
 src/openalea/hydroroot/model.py | 42 ++++++++++++++++-----------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/src/openalea/hydroroot/model.py b/src/openalea/hydroroot/model.py
index a8c881a..fdeb50a 100644
--- a/src/openalea/hydroroot/model.py
+++ b/src/openalea/hydroroot/model.py
@@ -22,62 +22,62 @@ class HydroRootModel(Model):
     length_data: list = declare(default=None, unit="[adim,m]", unit_comment="relative distance to the tip, meter",
                                 description="list of float or list of 2 list of float, lateral length vs relat. dist. to tip",
                                 min_value="", max_value="", value_comment="", references="", DOI="",
-                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     primary_length: float = declare(default=0.16, unit="m", unit_comment="", description="length of the primary root",
                                     min_value="0", max_value="", value_comment="", references="", DOI="",
-                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     seed: int = declare(default=None, unit="", unit_comment="", description="seed used to generate the architecture",
                         min_value="0", max_value="", value_comment="", references="", DOI="",
-                        variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                        variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     branching_delay: float = declare(default=2.0e-3, unit="m", unit_comment="",
                                      description="distance between branching laterals",
                                      min_value="0", max_value="", value_comment="", references="", DOI="",
-                                     variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                     variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     branching_variability: float = declare(default=0.25, unit="", unit_comment="",
                                            description="variability of the branching laterals distance",
                                            min_value="0", max_value="1", value_comment="", references="", DOI="",
-                                           variable_type="parameter", by="RsaModel", state_variable_type="",
+                                           variable_type="parameter", by="HydroRootModel", state_variable_type="",
                                            edit_by="user")
 
     order_max: int = declare(default=4, unit="", unit_comment="", description="maximum order of laterals",
                              min_value="0", max_value="1", value_comment="", references="", DOI="",
-                             variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                             variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
                                     min_value="0", max_value="", value_comment="", references="", DOI="",
-                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     nude_length: float = declare(default=0.02, unit="m", unit_comment="",
                                  description="part of roots without any lateral root, distance from tip",
                                  min_value="0", max_value="", value_comment="", references="", DOI="",
-                                 variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                 variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     ref_radius: float = declare(default=7.0e-5, unit="m", unit_comment="",
                                 description="reference radius or radius of the primary root",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
-                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     order_decrease_factor: float = declare(default=0.7, unit="", unit_comment="",
                                            description="radius decrease factor applied when increasing order",
                                            min_value="0", max_value="1", value_comment="", references="", DOI="",
-                                           variable_type="parameter", by="RsaModel", state_variable_type="",
+                                           variable_type="parameter", by="HydroRootModel", state_variable_type="",
                                            edit_by="user")
 
     # hydraulic parameters
     axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
-                                variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]",
                                             unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
-                                            variable_type="parameter", by="RsaModel", state_variable_type="",
+                                            variable_type="parameter", by="HydroRootModel", state_variable_type="",
                                             edit_by="user")
 
     # Soil1DModel input
@@ -151,11 +151,11 @@ class HydroRootModel(Model):
 
     surface: float = declare(default=None, unit="m2", unit_comment="", description="total surface of the RSA",
                              min_value="0", max_value="", value_comment="", references="", DOI="",
-                             variable_type="plant_scale_state", by="RsaModel", state_variable_type="", edit_by="user")
+                             variable_type="plant_scale_state", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     volume: float = declare(default=None, unit="m3", unit_comment="", description="total surface of the RSA",
                             min_value="0", max_value="", value_comment="", references="", DOI="",
-                            variable_type="plant_scale_state", by="RsaModel", state_variable_type="", edit_by="user")
+                            variable_type="plant_scale_state", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     def __init__(self, g, time_step, **scenario):
         # 4/21/26 FB: copied from root_cynaps RootWaterModel
@@ -359,20 +359,20 @@ def solute_solver_flux(self):
 
 @dataclass
 class HydrostaticModel(Model):
-    # Hydraulic input from RsaModel
+    # Hydraulic input from HydroRootModel
     length: float = declare(default=1.0e-4, unit="m", unit_comment="", description="vertices lentgh",
                                     min_value="0", max_value="", value_comment="", references="", DOI="",
-                                    variable_type="parameter", by="RsaModel", state_variable_type="", edit_by="user")
+                                    variable_type="parameter", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
-                                variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
+                                variable_type="input", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
-                                            variable_type="input", by="RsaModel", state_variable_type="",
+                                            variable_type="input", by="HydroRootModel", state_variable_type="",
                                             edit_by="user")
 
     # Soil1DModel input
@@ -423,16 +423,16 @@ def __init__(self, g, time_step, **scenario):
 
 @dataclass
 class SoluteModel(Model):
-    # Hydraulic input from RsaModel
+    # Hydraulic input from HydroRootModel
     axial_conductance_data: list = declare(default=None, unit="[m,10-9 m4.MPa-1.s-1]", unit_comment="distance to the tip, K",
                                 description="(2 list of Float) axial conductivity versus dist. to tip",
                                 min_value="0", max_value="", value_comment="", references="", DOI="",
-                                variable_type="input", by="RsaModel", state_variable_type="", edit_by="user")
+                                variable_type="input", by="HydroRootModel", state_variable_type="", edit_by="user")
 
     k0: list = declare(default=None, unit="[m,10-9 m.MPa-1.s-1]", unit_comment="distance to the tip, k",
                                             description="(2 list of Float) radial conductivity versus dist. to tip",
                                             min_value="0", max_value="", value_comment="", references="", DOI="",
-                                            variable_type="input", by="RsaModel", state_variable_type="",
+                                            variable_type="input", by="HydroRootModel", state_variable_type="",
                                             edit_by="user")
 
     # Soil1DModel input

From 0cc5429cb17585f7be164643304359cb5ed4e7e0 Mon Sep 17 00:00:00 2001
From: baugetfa <fabrice.bauget@cirad.fr>
Date: Mon, 27 Apr 2026 10:34:34 +0200
Subject: [PATCH 16/16] WIP to trigger CI

---
 doc/example/metafspm/example_metafspm_archi_from_file.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/doc/example/metafspm/example_metafspm_archi_from_file.py b/doc/example/metafspm/example_metafspm_archi_from_file.py
index d90e2d9..9216abe 100644
--- a/doc/example/metafspm/example_metafspm_archi_from_file.py
+++ b/doc/example/metafspm/example_metafspm_archi_from_file.py
@@ -160,9 +160,7 @@
 from openalea.hydroroot.init_parameter import Parameters
 parameter = Parameters()
 parameter.read_file('parameters_plant_01.yml')
-scenario = {}
-for key, value in parameter.__dict__.items():
-    scenario.update(value)
+scenario = parameter.metafspm_scenario()
 hydromodel = model.HydroRootModel(g, time_step, **scenario)
 hydromodel.hydrostatic_solver_flux()
 print('sap flux (microL/s): ', hydromodel.Jv)