diff --git a/pyposeidon/boundary.py b/pyposeidon/boundary.py
index f11f7e0..7ba181e 100644
--- a/pyposeidon/boundary.py
+++ b/pyposeidon/boundary.py
@@ -203,7 +203,7 @@ def buffer_(coasts, cbuffer):
         ww_.loc[idx] = b
 
     # join
-    wu = ww_.unary_union
+    wu = ww_.union_all()
     wu = gp.GeoDataFrame(geometry=[wu]).explode(index_parts=True).droplevel(0).reset_index(drop=True)
 
     rings = wu.boundary.is_ring  # get multiple boundaries instance
@@ -309,7 +309,7 @@ def tag(geometry, coasts, cbuffer, blevels):
             else:
                 try:
                     gw = gp.GeoDataFrame(
-                        geometry=list(ww.buffer(0).unary_union)
+                        geometry=list(ww.buffer(0).union_all())
                     )  # merge the polygons that are split (around -180/180)
                 except:
                     gw = gp.GeoDataFrame(geometry=list(ww.values))
@@ -342,7 +342,7 @@ def tag(geometry, coasts, cbuffer, blevels):
         block = coasts.copy()
 
     if not block.empty:
-        g = block.unary_union.symmetric_difference(grp)  # get the dif from the coasts
+        g = block.union_all().symmetric_difference(grp)  # get the dif from the coasts
     else:
         g = grp
 
@@ -529,9 +529,12 @@ def global_tag(geo, cbuffer, blevels, R=1):
 
     ww = gp.GeoDataFrame(geometry=cs)
 
-    gw = gp.GeoDataFrame(
-        geometry=list(ww.buffer(0).unary_union.geoms)
-    )  # merge the polygons that are split (around -180/180)
+    ww_ = ww.buffer(0).union_all()
+
+    if ww_.geom_type == "Polygon":
+        gw = gp.GeoDataFrame(geometry=[ww_])
+    else:
+        gw = gp.GeoDataFrame(geometry=list(ww_.geoms))  # merge the polygons that are split (around -180/180)
 
     gw = gp.GeoDataFrame(geometry=gw.boundary.values)
 
diff --git a/pyposeidon/dem_tools.py b/pyposeidon/dem_tools.py
index 66c6541..01de4ea 100644
--- a/pyposeidon/dem_tools.py
+++ b/pyposeidon/dem_tools.py
@@ -175,7 +175,7 @@ def fix(dem, coastline, **kwargs):
         block = shp.cx[minlon:maxlon, minlat:maxlat]
 
     try:
-        block = gp.GeoDataFrame(geometry=list(block.unary_union.geoms))
+        block = gp.GeoDataFrame(geometry=list(block.union_all().geoms))
     except:
         pass
 
@@ -189,7 +189,7 @@ def fix(dem, coastline, **kwargs):
     grp = grp.buffer(0.5)  # buffer it to get also the boundary points
 
     try:
-        g = block.unary_union.symmetric_difference(grp)  # get the diff
+        g = block.union_all().symmetric_difference(grp)  # get the diff
     except:
         g = grp  # no land
 
@@ -208,7 +208,7 @@ def fix(dem, coastline, **kwargs):
     except:
         b = shapely.geometry.GeometryCollection()
 
-    b = b.unary_union
+    b = b.union_all()
 
     # define wet/dry
     water = b
@@ -497,7 +497,7 @@ def check2(dataset, coastline):
     else:
         coasts = coastline
 
-    cc = coasts.unary_union
+    cc = coasts.union_all()
 
     wn = tree.query(cc, predicate="intersects").tolist()
 
diff --git a/pyposeidon/mgmsh.py b/pyposeidon/mgmsh.py
index 117d249..f74514c 100644
--- a/pyposeidon/mgmsh.py
+++ b/pyposeidon/mgmsh.py
@@ -166,8 +166,8 @@ def read_msh(filename, **kwargs):
         # nodes, tria = orient(nodes, tria, x="x", y="y")
         #        else:
         bgmesh = kwargs.get("bgmesh", None)
-        if not bgmesh:
-            tria = tria.reindex(columns=["a", "c", "b"])
+    #        if not bgmesh:
+    #            tria = tria.reindex(columns=["a", "c", "b"])
 
     # check if global and reproject
     if gglobal:  # convert to lat/lon
diff --git a/pyposeidon/utils/cast.py b/pyposeidon/utils/cast.py
index aedb385..d5f2252 100644
--- a/pyposeidon/utils/cast.py
+++ b/pyposeidon/utils/cast.py
@@ -280,7 +280,7 @@ def run(self, **kwargs):
 
         copy = get_value(self, kwargs, "copy", False)
 
-        ihot = get_value(self, kwargs, "ihot", 1)
+        ihot = get_value(self, kwargs, "ihot", 2)
 
         pwd = os.getcwd()
 
diff --git a/tests/test_schism_cast.py b/tests/test_schism_cast.py
index 2d054cc..8a1d0c4 100644
--- a/tests/test_schism_cast.py
+++ b/tests/test_schism_cast.py
@@ -198,7 +198,7 @@ def test_schism_cast_workflow(tmpdir):
         if not total.Dataset[var].equals(output.Dataset[var]):
             rb.append(var)
 
-    print(rb)
+    #    print(rb)
 
     #    flag = True TODO
     #    for var in rb:
@@ -208,5 +208,5 @@ def test_schism_cast_workflow(tmpdir):
     #            flag = True
     #    print(mdif)
 
-    expected = ["wetdry_side", "wetdry_elem", "wetdry_node", "zcor", "elev", "hvel"]
-    assert rb == expected
+    #    expected = ["wetdry_side", "wetdry_elem", "wetdry_node", "zcor", "elev", "hvel"]
+    assert (rb == ["zcor"]) or (rb == [])
diff --git a/tests/test_schism_reforecast.py b/tests/test_schism_reforecast.py
index 28d9605..1675a37 100644
--- a/tests/test_schism_reforecast.py
+++ b/tests/test_schism_reforecast.py
@@ -157,7 +157,7 @@ def test_schism_reforecast_workflow(tmpdir):
         if not total.Dataset[var].equals(r[var]):
             rb.append(var)
 
-    print(rb)
+    #    print(rb)
 
     #    flag = True TODO
     #    for var in rb:
@@ -166,5 +166,5 @@ def test_schism_reforecast_workflow(tmpdir):
     #        if mdif < 1.e-14 :
     #            flag = True
     #    print(mdif)
-    expected = ["wetdry_side", "wetdry_elem", "wetdry_node", "zcor", "elev", "hvel"]
-    assert rb == expected
+    #    expected = ["wetdry_side", "wetdry_elem", "wetdry_node", "zcor", "elev", "hvel"]
+    assert (rb == ["zcor"]) or (rb == [])