Use param.nml to define number of sediment classes

Before, it was hard-wired to 3.

Author: HansKimDWR

schimpy/schism_hotstart.py

@@ -306,6 +306,11 @@ def generate_3D_field(self, variable):
                 params = read_param_nml(self.sediment_fn)
                 self.Nbed = params['Nbed']
             kwargs_options.update({'Nbed': self.Nbed})
+
+            if 'nsed' not in self.__dict__.keys():
+                params = read_param_nml(self.param_nml)
+                self.nsed = params['sed_class']
+            kwargs_options.update({'nsed': self.nsed})
         if kwargs_options:
             var = VariableField(v_meta, variable, self.mesh,
                                 self.depths, self.date, self.crs,
@@ -468,8 +473,10 @@ def __init__(self, v_meta, vname, mesh, depths, date, crs, tr_mname, hotstart,
 
         if vname in ['SED3D_bed', 'SED3D_bedfrac']:
             self.Nbed = kwargs['Nbed']
+            self.nsed = kwargs['nsed']
         else:
             self.Nbed = 1
+            self.nsed = 1
         self.grid = self.get_grid()  # grid can be nodes/elems/edges
         self.n_hgrid = list(self.grid.values())[0][0]
         self.n_vgrid = list(self.grid.values())[1][0]
@@ -529,7 +536,7 @@ def get_grid(self):
                                      'MBEDP': (3, ['layer thickness', 'layer age', 'layer porosity'])},
                              'bedfrac': {'elem': (self.n_elems, self.elem),
                                          'Nbed': (self.Nbed, 0),
-                                         'nsed': (3, 0)}}
+                                         'nsed': (self.nsed, 0)}}
 
         if self.centering == 'edge':
             centering_options['edge'] = {'side': (self.n_edges, self.edge),
@@ -568,7 +575,7 @@ def define_new_grid(self, mesh):
                                      'MBEDP': (3, ['layer thickness', 'layer age', 'layer porosity'])},
                              'bedfrac': {'elem': (mesh.n_elems, mesh_elem),
                                          'Nbed': (self.Nbed, 0),
-                                         'nsed': (3, 0)}}
+                                         'nsed': (self.nsed, 0)}}
         if self.centering == 'edge':
             edge_depths = (mesh_depth[mesh.edges[:, 0]] +
                            mesh_depth[mesh.edges[:, 1]])/2
@@ -669,6 +676,13 @@ def simple_trend(self, ini_meta=None, inpoly=None):
         else:
             n_hgrid = self.n_hgrid
 
+        if self.variable_name in ['SED3D_bedfrac']:
+            value = str(value) # if a single float value is given, it needs to be converted to a string
+
+            # trigger an error if number of values is not equal to nsed
+            if len(value.split(',')) != self.nsed:
+                raise ValueError(
+                    "Number of values in %s does not match with nsed (%d)" % (value, self.nsed))
         if isinstance(value, (float, int)):
             # initialize the 3D field
             v_ini = np.full((n_hgrid, self.n_vgrid), value)