diff --git a/creating_geometric_models/Face.py b/creating_geometric_models/Face.py
index d6cf998..44e5ec3 100644
--- a/creating_geometric_models/Face.py
+++ b/creating_geometric_models/Face.py
@@ -1,15 +1,24 @@
 from collections import defaultdict
 import numpy as np
+from tqdm import tqdm, trange
+import pyvista as pv
 
 
 class Face:
     def __init__(self, vertex, connect):
-
         self.connect = connect
         self.vertex = vertex
         self.local_vid_lookup = {}
         self.ntriangles = self.connect.shape[0]
 
+    @classmethod
+    def from_pyvista(self, fname):
+        mesh = pv.read(fname)
+        mesh = mesh.clean(tolerance=0.1)
+        faces = mesh.faces.reshape((-1, 4))[:, 1:]
+        myFace = Face(mesh.points, faces)
+        return myFace
+
     @classmethod
     def from_ts(self, fid):
         surface_name = "undefined"
@@ -72,7 +81,9 @@ def proj(self, sProj):
 
         transformer = Transformer.from_crs("epsg:4326", sProj, always_xy=True)
         print("projecting the node coordinates")
-        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1])
+        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(
+            self.vertex[:, 0], self.vertex[:, 1]
+        )
         print(self.vertex)
         print("done projecting")
 
@@ -82,7 +93,9 @@ def convert_projected_to_latlon(self, sProj):
 
         transformer = Transformer.from_crs(sProj, "epsg:4326", always_xy=True)
         print("convert the node coordinates to lat/lon")
-        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(self.vertex[:, 0], self.vertex[:, 1])
+        self.vertex[:, 0], self.vertex[:, 1] = transformer.transform(
+            self.vertex[:, 0], self.vertex[:, 1]
+        )
         print(self.vertex)
         print("done converting")
 
@@ -109,11 +122,15 @@ def enforce_min_depth(self, depth_min):
         mesh = trimesh.Trimesh(self.vertex, self.connect)
         assert depth_min > 0
         # list vertex of the face boundary
-        unique_edges = mesh.edges[trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)]
-        if len(unique_edges)==0:
-            raise ValueError("the surface mesh has no boundary edges. "
+        unique_edges = mesh.edges[
+            trimesh.grouping.group_rows(mesh.edges_sorted, require_count=1)
+        ]
+        if len(unique_edges) == 0:
+            raise ValueError(
+                "the surface mesh has no boundary edges. "
                 "This means you extracted an incorrect set of surfaces "
-                "(the surface mesh of a closed region).")
+                "(the surface mesh of a closed region)."
+            )
         boundary_vid = np.array(list(set(list(unique_edges.flatten()))))
         # mask shallow water region
         shallow = np.where(mesh.vertices[:, 2] > -depth_min)
@@ -142,17 +159,38 @@ def __writeTs(self, fname, write_full_vertex_array=True, append=False):
 
         mode = "a" if append else "w"
         with open(fname, mode) as fout:
-            fout.write("GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n" % (fname))
-            for ivx in vertex_id_2_write:
-                fout.write("VRTX %s %s %s %s\n" % (ivx, self.vertex[ivx - 1, 0], self.vertex[ivx - 1, 1], self.vertex[ivx - 1, 2]))
-
-            for i in range(self.ntriangles):
-                fout.write("TRGL %d %d %d\n" % (self.connect[i, 0] + 1, self.connect[i, 1] + 1, self.connect[i, 2] + 1))
+            fout.write(
+                "GOCAD TSURF 1\nHEADER {\nname:%s\nborder: true\nmesh: false\n*border*bstone: true\n}\nTFACE\n"
+                % (fname)
+            )
+            for ivx in tqdm(vertex_id_2_write):
+                fout.write(
+                    "VRTX %s %s %s %s\n"
+                    % (
+                        ivx,
+                        self.vertex[ivx - 1, 0],
+                        self.vertex[ivx - 1, 1],
+                        self.vertex[ivx - 1, 2],
+                    )
+                )
+
+            for i in trange(self.ntriangles):
+                fout.write(
+                    "TRGL %d %d %d\n"
+                    % (
+                        self.connect[i, 0] + 1,
+                        self.connect[i, 1] + 1,
+                        self.connect[i, 2] + 1,
+                    )
+                )
             fout.write("END\n")
 
     def __computeNormal(self):
         "compute efficiently the normals"
-        normal = np.cross(self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :], self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :]).astype(float)
+        normal = np.cross(
+            self.vertex[self.connect[:, 1], :] - self.vertex[self.connect[:, 0], :],
+            self.vertex[self.connect[:, 2], :] - self.vertex[self.connect[:, 0], :],
+        ).astype(float)
         norm = np.linalg.norm(normal, axis=1)
         self.normal = np.divide(normal, norm[:, None])
 
@@ -161,11 +199,14 @@ def __writeStl(self, fname, append=False):
         mode = "a" if append else "w"
         with open(fname, mode) as fout:
             fout.write(f"solid {fname}\n")
-            for k in range(self.ntriangles):
+            for k in trange(self.ntriangles):
                 fout.write("facet normal %e %e %e\n" % tuple(self.normal[k, :]))
                 fout.write("outer loop\n")
                 for i in range(0, 3):
-                    fout.write("vertex %.10e %.10e %.10e\n" % tuple(self.vertex[self.connect[k, i], :]))
+                    fout.write(
+                        "vertex %.10e %.10e %.10e\n"
+                        % tuple(self.vertex[self.connect[k, i], :])
+                    )
                 fout.write("endloop\n")
                 fout.write("endfacet\n")
             fout.write(f"endsolid {fname}\n")
@@ -177,7 +218,7 @@ def __writebStl(self, fname):
         fout.seek(80)
         fout.write(struct.pack("<L", self.ntriangles))
         self.__computeNormal()
-        for k in range(self.ntriangles):
+        for k in trange(self.ntriangles):
             fout.write(struct.pack("<3f", *self.normal[k, :]))
             for i in range(0, 3):
                 fout.write(struct.pack("<3f", *self.vertex[self.connect[k, i], :]))
diff --git a/creating_geometric_models/Grid.py b/creating_geometric_models/Grid.py
index e3d998f..f0b3c2b 100644
--- a/creating_geometric_models/Grid.py
+++ b/creating_geometric_models/Grid.py
@@ -1,7 +1,7 @@
 import numpy as np
 import os
 from netCDF4 import Dataset
-
+from tqdm import tqdm
 
 class Grid:
     def __init__(self, fname, downsample):
@@ -53,9 +53,22 @@ def read_netcdf(self, fname, downsample):
         xvar = self.determine_netcdf_variables(["lon", "x"])
         yvar = self.determine_netcdf_variables(["lat", "y"])
         zvar = self.determine_netcdf_variables(["elevation", "z", "Band1"])
+
+        zvar_obj = fh.variables[zvar]
+        units = getattr(zvar_obj, "units", "").lower()
+
+        scale = 1.0  # default: meters
+
+        if units in ["m", "meter", "meters", "metre", "metres"]:
+            scale = 1.0
+        elif units in ["km", "kilometer", "kilometers"]:
+            scale = 1e3
+        else:
+            print(f"Warning: unknown elevation unit '{units}', assuming meters")
+
         self.x = fh.variables[xvar][0::downsample]
         self.y = fh.variables[yvar][0::downsample]
-        self.z = fh.variables[zvar][0::downsample, 0::downsample].astype(float)
+        self.z = fh.variables[zvar][0::downsample, 0::downsample].astype(float) * scale
 
         # transpose z if z was given as (lon, lat)
         dim_name_y = fh.variables[yvar].dimensions
@@ -64,7 +77,7 @@ def read_netcdf(self, fname, downsample):
             self.z = self.z.T
 
         if np.ma.is_masked(self.z):
-            self.z = np.ma.filled(self.z, float("nan")) * 1e3
+            self.z = np.ma.filled(self.z, float("nan"))
             self.is_sparse = True
         self.compute_nx_ny()
 
@@ -194,7 +207,7 @@ def isolate_hole(self, argHole):
         if argHole:
             x0, x1, y0, y1 = argHole
             print("tagging hole...")
-            for k in range(nconnect):
+            for k in tqdm(range(nconnect)):
                 coords = self.vertex[self.connect[k, :], 0:2]
                 xmin = min(coords[:, 0])
                 xmax = max(coords[:, 0])
@@ -208,12 +221,12 @@ def isolate_hole(self, argHole):
                 raise ValueError("no hole was tagged")
             print("done tagging hole")
 
-    def smooth(self, zrange):
+    def smooth(self, zrange, sig):
         "Smooth bathymetry in range +- zrange"
         "This is useful when preprocessing bathymetry data before intersection"
         import scipy.ndimage as ndimage
 
-        z_smooth = ndimage.gaussian_filter(self.z, sigma=(1, 1), order=0)
+        z_smooth = ndimage.gaussian_filter(self.z, sigma=(sig, sig), order=0)
         ids = np.where(abs(self.z) < zrange)
         self.z[ids] = z_smooth[ids]
 
diff --git a/creating_geometric_models/create_surface_from_rectilinear_grid.py b/creating_geometric_models/create_surface_from_rectilinear_grid.py
index 7fff090..916b731 100755
--- a/creating_geometric_models/create_surface_from_rectilinear_grid.py
+++ b/creating_geometric_models/create_surface_from_rectilinear_grid.py
@@ -9,8 +9,8 @@
 parser = argparse.ArgumentParser(description="create surface from a (possibly sparse, e.g. Slab2.0 dataset) structured dataset (e.g. netcdf)")
 parser.add_argument("input_file", help="netcdf file")
 parser.add_argument("output_file", help="output file (ext in stl, bstl, ts)")
-parser.add_argument("--subsample", nargs=1, type=int, metavar=("onesample_every"), default=[1], help="use only one value every onesample_every in both direction")
-parser.add_argument("--objectname", nargs=1, metavar=("objectname"), help="name of the surface in gocad")
+parser.add_argument("--subsample" type=int, metavar=("onesample_every"), default=[1], help="use only one value every onesample_every in both direction")
+parser.add_argument("--objectname", metavar=("objectname"), help="name of the surface in gocad")
 parser.add_argument("--hole", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="isolate a hole in surface defined by x0<=x<=x1 and y0<=y<=y1 (stl and ts output only)", type=float)
 parser.add_argument("--crop", nargs=4, metavar=(("x0"), ("x1"), ("y0"), ("y1")), help="select only surfaces in x0<=x<=x1 and y0<=y<=y1", type=float)
 parser.add_argument(
@@ -21,23 +21,23 @@
  projname: name of the (projected) Coordinate Reference System (CRS) (e.g. EPSG:32646 for UTM46N)",
 )
 parser.add_argument("--translate", nargs=2, metavar=("x0", "y0"), default=([0, 0]), help="translates all nodes by (x0,y0)", type=float)
-parser.add_argument("--smooth", nargs=1, metavar=("zrange"), help="smooth zgrid for -zrange < z < zrange", type=float)
-parser.add_argument("--change_zero_elevation", nargs=1, metavar=("new_elevation"), help="change elevation of points with z=0, thus avoiding errors when interesecting with sea surface", type=float)
+parser.add_argument("--smooth", nargs=2, metavar=("zrange", "sigma"), help="smooth zgrid for -zrange < z < zrange", type=float)
+parser.add_argument("--change_zero_elevation", metavar=("new_elevation"), help="change elevation of points with z=0, thus avoiding errors when interesecting with sea surface", type=float)
 args = parser.parse_args()
 
 if not args.objectname:
     base = os.path.basename(args.input_file)
     args.objectname = os.path.splitext(base)[0]
 else:
-    args.objectname = args.objectname[0]
+    args.objectname = args.objectname
 
-structured_grid = Grid(args.input_file, args.subsample[0])
+structured_grid = Grid(args.input_file, args.subsample)
 structured_grid.crop(args.crop)
 
 if args.change_zero_elevation:
-    structured_grid.change_zero_elevation(args.change_zero_elevation[0])
+    structured_grid.change_zero_elevation(args.change_zero_elevation)
 if args.smooth:
-    structured_grid.smooth(args.smooth[0])
+    structured_grid.smooth(args.smooth[0], rgs.smooth[1])
 
 structured_grid.generate_vertex()
 
@@ -62,7 +62,7 @@
     for sid in range(nsolid):
         idtr = np.where(structured_grid.solid_id == sid)[0]
         aVid = np.unique(structured_grid.connect[idtr, :].flatten())
-        myFace = Face(vertex=None, connect=structured_grid.connect[idtr, :])
+        myFace = Face(vertex=structured_grid.vertex, connect=structured_grid.connect[idtr, :])
         if structured_grid.is_sparse:
             myFace.reindex(structured_grid.vid_lookup)
-        myFace.write(f"{basename}{sid}{ext}", structured_grid.vertex, write_full_vertex_array=False)
+        myFace.write(f"{basename}{sid}{ext}", write_full_vertex_array=False)
diff --git a/creating_geometric_models/refine_and_smooth_mesh.py b/creating_geometric_models/refine_and_smooth_mesh.py
index 635ad31..a935e4b 100755
--- a/creating_geometric_models/refine_and_smooth_mesh.py
+++ b/creating_geometric_models/refine_and_smooth_mesh.py
@@ -3,6 +3,7 @@
 import argparse
 import os
 import trimesh
+import numpy as np
 
 parser = argparse.ArgumentParser(description="refine and smooth")
 parser.add_argument("input_file", help="surface in ts file format")
@@ -22,19 +23,81 @@
     required=True,
     help="number of and refine/smoothing steps",
 )
+
+parser.add_argument(
+    "--remesh_first",
+    nargs=1,
+    metavar=("mesh_size"),
+    type=float,
+    help="remesh first the surface with pygalmesh",
+)
+
+parser.add_argument(
+    "--fix_boundary",
+    dest="fix_boundary",
+    action="store_true",
+    help="fix boundary edge when interpolating",
+)
+
+
 args = parser.parse_args()
 
 fid = open(args.input_file)
 myFace = Face.from_ts(fid)
 
-a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect)
+if args.remesh_first:
+    """
+    # pygalmesh installed with:
+    sudo apt install libcgal-dev libeigen3-dev
+    pip install pygalmesh
+    """
+    import pygalmesh
+    import meshio
+
+    mesh = meshio.Mesh(points=myFace.vertex, cells=[("triangle", myFace.connect)])
+    mesh.write("trash_me.vtk")
+    mesh_size = args.remesh_first[0]
+    print(f"remeshing with pygalmesh aiming for mesh size {mesh_size}")
+    mesh = pygalmesh.remesh_surface(
+        "trash_me.vtk",
+        max_edge_size_at_feature_edges=mesh_size,
+        min_facet_angle=25,
+        max_radius_surface_delaunay_ball=mesh_size,
+        max_facet_distance=mesh_size,
+        verbose=False,
+    )
+    os.remove("trash_me.vtk")
+    print("done remeshing")
+    a = trimesh.Trimesh(vertices=mesh.points, faces=mesh.cells[0].data)
+else:
+    a = trimesh.Trimesh(vertices=myFace.vertex, faces=myFace.connect)
+
 
 for i in range(args.N[0]):
     a = a.subdivide()
 for i in range(args.P[0]):
     a = a.subdivide()
+
+    if args.fix_boundary:
+        # Identify boundary edges using grouping
+        unique_edges = a.edges[
+            trimesh.grouping.group_rows(a.edges_sorted, require_count=1)
+        ]
+        # Extract unique vertices on boundary
+        boundary_vertices = np.unique(unique_edges)
+        a_before_smoothing = a.copy()
+
     a = trimesh.smoothing.filter_taubin(a)
 
+    if args.fix_boundary:
+        # Replace only the internal vertices
+        a.vertices[np.isin(np.arange(len(a.vertices)), boundary_vertices)] = (
+            a_before_smoothing.vertices[
+                np.isin(np.arange(len(a.vertices)), boundary_vertices)
+            ]
+        )
+
+
 myFace = Face(a.vertices, a.faces)
 basename, ext = os.path.splitext(args.input_file)
 myFace.write(f"{basename}_refined_smooth_{args.N[0]}{ext}")