Quhall error fix

Author: browniea

mhkit/river/io/d3d.py

@@ -650,6 +650,34 @@ def variable_interpolation(
 
     if not isinstance(to_pandas, bool):
         raise TypeError(f"to_pandas must be of type bool. Got: {type(to_pandas)}")
+    
+    #Avoid Quhall errors by filling in bad sigma values with -1 for the surface or 0 for the waterbed based on neighboring values
+    if "mesh2d_interface_sigma" in data.variables:
+        sigma = data["mesh2d_interface_sigma"].values
+        bad = ~np.isfinite(sigma) | (np.abs(sigma) > 0.1)
+        if bad.any():
+            # Check neighboring values to determine if sigma should be -1 or 0
+            for idx in np.where(bad)[0]:
+                neighbors = []
+                if idx > 0:
+                    neighbors.append(sigma[idx - 1])
+                if idx < len(sigma) - 1:
+                    neighbors.append(sigma[idx + 1])
+                
+                if neighbors:
+                    avg_neighbor = np.mean(neighbors)
+                    if np.abs(avg_neighbor - (-1.0)) < np.abs(avg_neighbor):
+                        sigma[idx] = -1.0
+                    else:
+                        sigma[idx] = 0.0
+                else:
+                    sigma[idx] = -1.0
+            
+            data["mesh2d_interface_sigma"] = xr.DataArray(
+                sigma,
+                dims=data["mesh2d_interface_sigma"].dims,
+                attrs=data["mesh2d_interface_sigma"].attrs,
+                )
 
     data_raw = {}
     for var in variables: