Zero layer thickness layers are considered active due to float roundi…

…ng (#324)

* Zero layer thickness layers are considered active due to float rounding

Subtracting two botms to calculate the thickness results is some layers with a thickness of ~1e-9. This rounding later results in a problem because this positive thickness leads to an idomain of 1.

This can either be solved in `nlmod.layers.calculating_thickness()` or in `nlmod.layers.get_idomain()`. What do you prefer?

* Remove trailing whitespace

* Applied patch also to get_idomain

nlmod/dims/layers.py

@@ -47,6 +47,10 @@ def calculate_thickness(ds, top="top", bot="botm"):
                     thickness[lay] = ds[bot][lay - 1] - ds[bot][lay]
         else:
             raise ValueError("2d top should have same last dimension as bot")
+
+    # subtracting floats can result in rounding errors. Mainly anoying for zero thickness layers.
+    thickness = thickness.where(~np.isclose(thickness, 0.), 0.)
+
     if isinstance(ds[bot], xr.DataArray):
         thickness.name = "thickness"
         if hasattr(ds[bot], "long_name"):
@@ -1310,6 +1314,8 @@ def get_idomain(ds):
     idomain.attrs.clear()
     # set idomain of cells  with a positive thickness to 1
     thickness = calculate_thickness(ds)
+    # subtracting floats can result in rounding errors. Mainly anoying for zero thickness layers.
+    thickness = thickness.where(~np.isclose(thickness, 0.), 0.)
     idomain.data[thickness.data > 0.0] = 1
     # set idomain above/below the first/last active layer to 0
     idomain.data[idomain.where(idomain > 0).ffill(dim="layer").isnull()] = 0