RubarBE: split critical Shields stress (TCADIM) from water viscosity (VISC)

src/Solver/RubarBE.py

@@ -92,6 +92,7 @@ class Rubar3(CommandLineSolver):
             ("rubarbe_sediment_mult_3", ""),
             ("rubarbe_sediment_mult_4", ""),
             ("rubarbe_sediment_mult_5", ""),
+            ("rubarbe_sediment_tcadim", "0.047"),
             ("rubarbe_sediment_visc", "1.3E-06"),
             ("rubarbe_sediment_opts", "3"),
             ("rubarbe_sediment_odchar", "0"),
@@ -200,7 +201,7 @@ class Rubar3(CommandLineSolver):
     _ts_fields = (
         ("ros", "real"), ("por", "real"), ("dcharg", "real"),
         ("halfa", "real"), ("mult_1", "multipliers"),
-        ("visc", "real"), ("opts", "law"), ("odchar", "integer"),
+        ("visc", "visc"), ("opts", "law"), ("odchar", "integer"),
         ("unisol", "integer"), ("typdef", "integer"),
         ("depot", "depot"), ("choixc", "integer"),
         ("option", "options"), ("capsol", "integer"),
@@ -208,6 +209,28 @@ class Rubar3(CommandLineSolver):
         ("defond", "integer"), ("varcons", "integer"),
     )
 
+    # rubar_init.f90 (inits): laws whose critical stress may come from the
+    # Shields diagram, and laws without critical stress
+    _shields_laws = {1, 10, 11, 12, 14, 25, 26}
+    _no_threshold_laws = {2, 3, 8, 13, 18, 19, 22, 23, 24}
+    _brownlie_law = 16
+
+    @classmethod
+    def tcadim_required(cls, opts, option):
+        """True if RubarBE reads the dimensionless critical stress TCADIM
+        for transport law OPTS and options OPTION (notice §3.10)"""
+        try:
+            law = int(opts)
+        except ValueError:
+            return False
+
+        shields = len(option) == 3 and option[1] != "1"
+        if shields and law in cls._shields_laws:
+            return False
+        if law in cls._no_threshold_laws or law == cls._brownlie_law:
+            return False
+        return True
+
     @staticmethod
     def _reaches_without_solid_boundary(study):
         river = study.river
@@ -344,6 +367,26 @@ class Rubar3(CommandLineSolver):
             value = self._required(values, key)
             if kind == "real":
                 result = self._real(value, key)
+            elif kind == "visc":
+                # Notice annex 1: TCADIM then VISC, read as (A50, A20) then
+                # two F10.0. A single value would be read as TCADIM, with
+                # VISC = 0 (rubar_init.f90, inits).
+                if Decimal(self._real(value, key)) <= 0:
+                    raise ValueError(f"{key}: expected a positive viscosity")
+                tcadim_key = prefix + "tcadim"
+                tcadim = values.get(tcadim_key, "")
+                if tcadim:
+                    result = self._real(tcadim, tcadim_key)
+                elif self.tcadim_required(values.get(prefix + "opts", ""),
+                                          values.get(prefix + "option", "")):
+                    message = trad["missing_tcadim"].format(key=tcadim_key)
+                    logger.warning(message)
+                    if qlog is not None:
+                        qlog.put(trad["warning"].format(message=message))
+                    return False
+                else:
+                    result = " " * 10
+                result += self._real(value, key)
             elif kind == "multipliers":
                 result = self._real(value, key)
                 extras = [values.get(prefix + f"mult_{i}", "")