diff --git a/test/grid_geoflow.exo b/test/grid_geoflow.exo
deleted file mode 100644
index 8abc65564..000000000
Binary files a/test/grid_geoflow.exo and /dev/null differ
diff --git a/test/test_fesom.py b/test/test_fesom.py
new file mode 100644
index 000000000..11f479834
--- /dev/null
+++ b/test/test_fesom.py
@@ -0,0 +1,13 @@
+import uxarray as ux
+
+import os
+from pathlib import Path
+
+current_path = Path(os.path.dirname(os.path.realpath(__file__)))
+
+fesom_ugrid_diag_file= current_path / "meshfiles" / "ugrid" / "fesom" / "fesom.mesh.diag.nc"
+
+
+def test_open_fesom_ugrid():
+    uxgrid = ux.open_grid(fesom_ugrid_diag_file)
+    uxgrid.validate()
diff --git a/uxarray/grid/connectivity.py b/uxarray/grid/connectivity.py
index 1fe9efbaf..ee40c45e6 100644
--- a/uxarray/grid/connectivity.py
+++ b/uxarray/grid/connectivity.py
@@ -58,8 +58,8 @@ def close_face_nodes(face_node_connectivity, n_face, n_max_face_nodes):
 
 
 def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
-    """Replaces all instances of the the current fill value (``original_fill``)
-    in (``grid_var``) with (``new_fill``) and converts to the dtype defined by
+    """Replaces all instances of the current fill value (``original_fill``) in
+    (``grid_var``) with (``new_fill``) and converts to the dtype defined by
     (``new_dtype``)
 
     Parameters
@@ -82,6 +82,7 @@ def _replace_fill_values(grid_var, original_fill, new_fill, new_dtype=None):
     # locations of fill values
     if original_fill is not None and np.isnan(original_fill):
         fill_val_idx = np.isnan(grid_var)
+        grid_var[fill_val_idx] = 0.0  # todo?
     else:
         fill_val_idx = grid_var == original_fill
 
diff --git a/uxarray/io/_ugrid.py b/uxarray/io/_ugrid.py
index 4af5959ef..3629a843b 100644
--- a/uxarray/io/_ugrid.py
+++ b/uxarray/io/_ugrid.py
@@ -11,12 +11,10 @@ def _read_ugrid(ds):
     """Parses an unstructured grid dataset and encodes it in the UGRID
     conventions."""
 
-    # Grid Topology
+    # Extract grid topology attributes
     grid_topology_name = list(ds.filter_by_attrs(cf_role="mesh_topology").keys())[0]
     ds = ds.rename({grid_topology_name: "grid_topology"})
 
-    # Coordinates
-
     # get the names of node_lon and node_lat
     node_lon_name, node_lat_name = ds["grid_topology"].node_coordinates.split()
     coord_dict = {
@@ -37,7 +35,6 @@ def _read_ugrid(ds):
         coord_dict[face_lat_name] = ugrid.FACE_COORDINATES[1]
 
     ds = ds.rename(coord_dict)
-    # Connectivity
 
     conn_dict = {}
     for conn_name in ugrid.CONNECTIVITY_NAMES:
@@ -56,25 +53,38 @@ def _read_ugrid(ds):
     dim_dict = {}
 
     # Rename Core Dims (node, edge, face)
-    if "node_dimension" in ds["grid_topology"]:
+    if "node_dimension" in ds["grid_topology"].attrs:
         dim_dict[ds["grid_topology"].node_dimension] = ugrid.NODE_DIM
     else:
         dim_dict[ds["node_lon"].dims[0]] = ugrid.NODE_DIM
 
-    if "face_dimension" in ds["grid_topology"]:
+    if "face_dimension" in ds["grid_topology"].attrs:
         dim_dict[ds["grid_topology"].face_dimension] = ugrid.FACE_DIM
     else:
         dim_dict[ds["face_node_connectivity"].dims[0]] = ugrid.FACE_DIM
 
-    if "edge_dimension" in ds["grid_topology"]:
+    if "edge_dimension" in ds["grid_topology"].attrs:
         # edge dimension is not always provided
         dim_dict[ds["grid_topology"].edge_dimension] = ugrid.EDGE_DIM
     else:
         if "edge_lon" in ds:
             dim_dict[ds["edge_lon"].dims[0]] = ugrid.EDGE_DIM
 
+    for conn_name in conn_dict.values():
+        # Ensure grid dimension (i.e. 'n_face') is always the first dimension
+        da = ds[conn_name]
+        dims = da.dims
+
+        for grid_dim in dim_dict.keys():
+            if dims[1] == grid_dim:
+                ds[conn_name] = da.T
+
     dim_dict[ds["face_node_connectivity"].dims[1]] = ugrid.N_MAX_FACE_NODES_DIM
 
+    for dim in ds.dims:
+        if ds.sizes[dim] == 2:
+            dim_dict[dim] = "two"
+
     ds = ds.swap_dims(dim_dict)
 
     return ds, dim_dict
@@ -149,7 +159,9 @@ def _standardize_connectivity(ds, conn_name):
         if "start_index" in ds[conn_name].attrs:
             new_conn[new_conn != INT_FILL_VALUE] -= INT_DTYPE(ds[conn_name].start_index)
         else:
-            new_conn[new_conn != INT_FILL_VALUE] -= INT_DTYPE(new_conn.min())
+            fill_value_indices = new_conn != INT_FILL_VALUE
+            start_index = new_conn[fill_value_indices].min()
+            new_conn[fill_value_indices] -= INT_DTYPE(start_index)
 
         # reassign data to use updated connectivity
         ds[conn_name].data = new_conn