diff --git a/flopy/mf6/data/mfdataarray.py b/flopy/mf6/data/mfdataarray.py
index 30752aed5..beeff7657 100644
--- a/flopy/mf6/data/mfdataarray.py
+++ b/flopy/mf6/data/mfdataarray.py
@@ -1602,10 +1602,10 @@ def _set_storage_netcdf(self, nc_dataset, create=False):
         input_tag = f"{m}/{p}/{n}".lower()
 
         if create:
-            # cache data and update file before resetting storage
+            # add array to netcdf dataset
             nc_varname = f"{p}_{n}".lower()
             data = self._get_data()
-            nc_dataset.create_var(nc_varname, input_tag, 0, data, self.structure.shape)
+            nc_dataset.create_array(nc_varname, input_tag, 0, data, self.structure.shape)
 
         storage._set_storage_netcdf(
             nc_dataset, input_tag, self.structure.layered, storage.layer_storage.get_total_size() - 1
diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py
index 075abe4fc..58ecb2084 100644
--- a/flopy/mf6/data/mffileaccess.py
+++ b/flopy/mf6/data/mffileaccess.py
@@ -473,7 +473,7 @@ def load_netcdf_array(
         nc_tag,
         layer,
     ):
-        return nc_dataset.data(nc_tag, layer + 1)
+        return nc_dataset.array(nc_tag, layer + 1)
 
     def get_data_string(self, data, data_type, data_indent=""):
         layer_data_string = [str(data_indent)]
diff --git a/flopy/mf6/data/mfstructure.py b/flopy/mf6/data/mfstructure.py
index e1ff3609b..7b127b1f8 100644
--- a/flopy/mf6/data/mfstructure.py
+++ b/flopy/mf6/data/mfstructure.py
@@ -1425,7 +1425,7 @@ def __init__(self, data_item, model_data, package_type, dfn_list):
             or "nodes" in data_item.shape
             or len(data_item.layer_dims) > 1
         )
-        # TODO: only gwf, gwt, gwe models
+        # TODO: only gwf, gwt, gwe models NETCDF-DEV
         self.netcdf = (
             ("ncol" in data_item.shape
             or "nrow" in data_item.shape
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index 033364da7..87d170adb 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -940,19 +940,26 @@ def load_base(
                 dis_str = {
                     "dis6": "structured",
                     "disv6": "vertex",
-                    "disu6": "unstructured",
+                    #"disu6": "unstructured",
                 }
                 dis_type = None
                 for t in instance.name_file.packages.get_data():
                     if t[0].lower().startswith("dis"):
                         dis_type = t[0].lower()
                         break
-                if dis_type:
+                if dis_type and dis_type in dis_str:
                     nc_fpth = os.path.join(instance.model_ws, nc_filerecord[0][0])
                     instance._nc_dataset = open_netcdf_dataset(nc_fpth, dis_type=dis_str[dis_type])
                 else:
-                    pass
-                    # TODO: set error "invalid dis for netcdf input"
+                    message = (
+                        "Invalid discretization type "
+                        f"provided for model {modelname} "
+                        "NetCDF input"
+                    )
+                    raise MFDataException(
+                        model=modelname,
+                        message=message,
+                    )
 
         # order packages
         vnum = mfstructure.MFStructure().get_version_string()
diff --git a/flopy/utils/model_netcdf.py b/flopy/utils/model_netcdf.py
index 3753826e7..4b9193504 100644
--- a/flopy/utils/model_netcdf.py
+++ b/flopy/utils/model_netcdf.py
@@ -46,11 +46,8 @@ def open(self, nc_fpth: str) -> None:
         self._dpth = fpth.parent
         self._name = fpth.name
 
-        try:
-            self._dataset = xr.open_dataset(fpth, engine="netcdf4")
-            self._set_mapping()
-        except Exception as e:
-            print(f"Exception: {e}")
+        self._dataset = xr.open_dataset(fpth, engine="netcdf4")
+        self._set_mapping()
 
     def create(
         self, model_type: str, model_name: str, nc_type: str, fname: str, dis: Grid
@@ -60,24 +57,20 @@ def create(
         self._name = fname
         self._tags = {}
 
-        assert self._type == "mesh2d" or self._type == "structured"
+        if self._type != "mesh2d" and self._type != "structured":
+            raise Exception('Supported NetCDF file types are "mesh2d" and "structured"')
         if isinstance(dis, VertexGrid) and self._type != "mesh2d":
-            # TODO error
-            pass
-
-        try:
-            self._dataset = xr.Dataset()
-            self._set_global_attrs(model_type, model_name)
-            self._set_grid(dis)
-        except Exception as e:
-            print(f"Exception: {e}")
+            raise Exception("VertexGrid object must use mesh2d netcdf file type")
 
+        self._dataset = xr.Dataset()
+        self._set_global_attrs(model_type, model_name)
+        self._set_grid(dis)
         # print(self._dataset.info())
 
-    def create_var(
+    def create_array(
         self, varname: str, tag: str, layer: int, data: np.typing.ArrayLike, shape: list
     ):
-        raise NotImplementedError("create_var not implemented in base class")
+        raise NotImplementedError("create_array not implemented in base class")
 
     def write(self, path: str) -> None:
         nc_fpath = Path(path) / self._name
@@ -90,14 +83,13 @@ def layered(self) -> bool:
 
         return res
 
-    def data(self, path, layer=None):
-        raise NotImplementedError("data not implemented in base class")
+    def array(self, path, layer=None):
+        raise NotImplementedError("array not implemented in base class")
 
     def _set_mapping(self):
         var_d = {}
         if ("modflow6_grid" and "modflow6_model") not in self._dataset.attrs:
-            # TODO error
-            pass
+            raise Exception("Invalid MODFLOW 6 netcdf dataset")
         else:
             self._modelname = (
                 self._dataset.attrs["modflow6_model"].split(":")[0].lower()
@@ -133,8 +125,7 @@ def _set_global_attrs(self, model_type, model_name):
             dep_var = "temperature"
             model = "Groundwater Energy (GWE)"
         else:
-            # TODO error?
-            pass
+            raise Exception("NetCDF supported for GWF, GWT and GWE models")
 
         if self._type == "structured":
             grid = self._type.upper()
@@ -155,48 +146,43 @@ def _set_global_attrs(self, model_type, model_name):
     def _set_grid(self, dis):
         raise NotImplementedError("_set_grid not implemented in base class")
 
-    def _create_var(
+    def _create_array(
         self, varname: str, tag: str, data: np.typing.ArrayLike, nc_shape: list
     ):
-        try:
-            layer = 0
-            var_d = {varname: (nc_shape, data)}
+        layer = 0
+        var_d = {varname: (nc_shape, data)}
+        self._dataset = self._dataset.assign(var_d)
+        self._dataset[varname].attrs["modflow6_input"] = tag
+        if tag not in self._tags:
+            self._tags[tag] = {}
+        if layer in self._tags[tag]:
+            raise Exception(f"Array variable tag already exists: {tag}")
+        self._tags[tag][layer] = varname
+
+    def _create_layered_array(
+        self, varname: str, tag: str, data: np.typing.ArrayLike, nc_shape: list
+    ):
+        for layer in range(data.shape[0]):
+            mf6_layer = layer + 1
+            layer_vname = f"{varname}_l{mf6_layer}"
+            var_d = {layer_vname: (nc_shape, data[layer].flatten())}
             self._dataset = self._dataset.assign(var_d)
-            self._dataset[varname].attrs["modflow6_input"] = tag
+            self._dataset[layer_vname].attrs["modflow6_input"] = tag
+            self._dataset[layer_vname].attrs["modflow6_layer"] = layer + 1
             if tag not in self._tags:
                 self._tags[tag] = {}
-            if layer in self._tags[tag]:
-                # TODO error?
-                pass
-            self._tags[tag][layer] = varname
-        except Exception as e:
-            print(f"Exception: {e}")
-
-    def _create_layered_var(
-        self, varname: str, tag: str, data: np.typing.ArrayLike, nc_shape: list
-    ):
-        try:
-            for layer in range(data.shape[0]):
-                layer_vname = f"{varname}_l{layer + 1}"
-                var_d = {layer_vname: (nc_shape, data[layer].flatten())}
-                self._dataset = self._dataset.assign(var_d)
-                self._dataset[layer_vname].attrs["modflow6_input"] = tag
-                self._dataset[layer_vname].attrs["modflow6_layer"] = layer + 1
-                if tag not in self._tags:
-                    self._tags[tag] = {}
-                if layer in self._tags[tag]:
-                    # TODO error?
-                    pass
-                self._tags[tag][layer + 1] = layer_vname
-        except Exception as e:
-            print(f"Exception: {e}")
+            if mf6_layer in self._tags[tag]:
+                raise Exception(
+                    f"Array variable tag already exists: {tag}, layer={layer}"
+                )
+            self._tags[tag][mf6_layer] = layer_vname
 
 
 class DisNetCDFStructured(ModelNetCDFDataset):
     def __init__(self):
         super().__init__()
 
-    def create_var(
+    def create_array(
         self, varname: str, tag: str, layer: int, data: np.typing.ArrayLike, shape: list
     ):
         nc_shape = None
@@ -209,25 +195,21 @@ def create_var(
                 nc_shape = ["y"]
             elif shape[0].lower() == "ncol":
                 nc_shape = ["x"]
-            else:
-                # TODO error?
-                pass
 
-        self._create_var(varname, tag, data, nc_shape)
+        self._create_array(varname, tag, data, nc_shape)
 
-    def data(self, path, layer=None):
+    def array(self, path, layer=None):
         if len(self._dataset[self._tags[path][0]].data.shape) == 3:
             return self._dataset[self._tags[path][0]].data[layer - 1]
         else:
             return self._dataset[self._tags[path][0]].data
 
     def _set_grid(self, dis):
-        x = dis.xcellcenters[0]
-        y = dis.ycellcenters[:, 0]
-        z = list(range(1, dis.nlay + 1))
+        # lenunits = {0: "undefined", 1: "feet", 2: "meters", 3: "centimeters"}
+        lenunits = {0: "m", 1: "feet", 2: "m", 3: "cm"}
 
         x_bnds = []
-        for idx, x in enumerate(dis.xyedges[0]):
+        for idx, val in enumerate(dis.xyedges[0]):
             if idx + 1 < len(dis.xyedges[0]):
                 bnd = []
                 bnd.append(dis.xyedges[0][idx])
@@ -235,13 +217,17 @@ def _set_grid(self, dis):
                 x_bnds.append(bnd)
 
         y_bnds = []
-        for idx, y in enumerate(dis.xyedges[1]):
+        for idx, val in enumerate(dis.xyedges[1]):
             if idx + 1 < len(dis.xyedges[1]):
                 bnd = []
                 bnd.append(dis.xyedges[1][idx + 1])
                 bnd.append(dis.xyedges[1][idx])
                 y_bnds.append(bnd)
 
+        x = dis.xcellcenters[0]
+        y = dis.ycellcenters[:, 0]
+        z = list(range(1, dis.nlay + 1))
+
         # create coordinate vars
         var_d = {"z": (["z"], z), "y": (["y"], y), "x": (["x"], x)}
         self._dataset = self._dataset.assign(var_d)
@@ -253,29 +239,25 @@ def _set_grid(self, dis):
         # set coordinate variable attributes
         self._dataset["z"].attrs["units"] = "layer"
         self._dataset["z"].attrs["long_name"] = "layer number"
-        self._dataset["y"].attrs["units"] = "m"  # TODO in dis?
+        self._dataset["y"].attrs["units"] = lenunits[dis.lenuni]
         self._dataset["y"].attrs["axis"] = "Y"
         self._dataset["y"].attrs["standard_name"] = "projection_y_coordinate"
         self._dataset["y"].attrs["long_name"] = "Northing"
-        self._dataset["y"].attrs["grid_mapping"] = "projection"  # TODO
+        # self._dataset["y"].attrs["grid_mapping"] = "projection"  # TODO
         self._dataset["y"].attrs["bounds"] = "y_bnds"
-        self._dataset["x"].attrs["units"] = "m"  # TODO in dis?
+        self._dataset["x"].attrs["units"] = lenunits[dis.lenuni]
         self._dataset["x"].attrs["axis"] = "X"
         self._dataset["x"].attrs["standard_name"] = "projection_x_coordinate"
         self._dataset["x"].attrs["long_name"] = "Easting"
-        self._dataset["x"].attrs["grid_mapping"] = "projection"  # TODO
+        # self._dataset["x"].attrs["grid_mapping"] = "projection"  # TODO
         self._dataset["x"].attrs["bounds"] = "x_bnds"
 
 
-# TODO: base class ModelNetCDFMesh2d(ModelNetCDFDataset) ?
-#       for common routines e.g. data()
-
-
 class DisNetCDFMesh2d(ModelNetCDFDataset):
     def __init__(self):
         super().__init__()
 
-    def create_var(
+    def create_array(
         self, varname: str, tag: str, layer: int, data: np.typing.ArrayLike, shape: list
     ):
         nc_shape = None
@@ -288,11 +270,11 @@ def create_var(
             nc_shape = ["nmesh_face"]
 
         if len(data.shape) == 3:
-            self._create_layered_var(varname, tag, data, nc_shape)
+            self._create_layered_array(varname, tag, data, nc_shape)
         else:
-            self._create_var(varname, tag, data.flatten(), nc_shape)
+            self._create_array(varname, tag, data.flatten(), nc_shape)
 
-    def data(self, path, layer=None):
+    def array(self, path, layer=None):
         if not layer:
             layer = 0
         if path in self._tags:
@@ -305,7 +287,7 @@ def _set_grid(self, dis):
         print(dir(dis))
         print(dis.ncpl)
         print(dis.nvert)
-        print(dis.get_cell_vertices())
+        # print(dis.get_cell_vertices())
         # nmesh_node = dis.nvert
         # nmesh_face = dis.ncpl
         # max_nmesh_face_nodes = 4 ; # assume 4 for dis?
@@ -316,17 +298,17 @@ class DisvNetCDFMesh2d(ModelNetCDFDataset):
     def __init__(self):
         super().__init__()
 
-    def create_var(
+    def create_array(
         self, varname: str, tag: str, layer: int, data: np.typing.ArrayLike, shape: list
     ):
         nc_shape = ["nmesh_face"]
 
         if len(data.shape) == 2:
-            self._create_layered_var(varname, tag, data, nc_shape)
+            self._create_layered_array(varname, tag, data, nc_shape)
         else:
-            self._create_var(varname, tag, data.flatten(), nc_shape)
+            self._create_array(varname, tag, data.flatten(), nc_shape)
 
-    def data(self, path, layer=None):
+    def array(self, path, layer=None):
         if not layer:
             layer = 0
         if path in self._tags:
@@ -343,35 +325,37 @@ def open_netcdf_dataset(nc_fpth: str, dis_type: str) -> ModelNetCDFDataset:
     # dis_type corresponds to a flopy.discretization derived object type
     nc_dataset = None
     dis_str = dis_type.lower()
-    assert dis_str == "structured" or dis_str == "vertex"
+    if dis_str != "vertex" and dis_str != "structured":
+        raise Exception(
+            "Supported NetCDF discretication types "
+            'are "vertex" (DISV) and "structured" '
+            "(DIS)"
+        )
 
     fpth = Path(nc_fpth).resolve()
-    try:
-        dataset = xr.open_dataset(fpth, engine="netcdf4")
-        if ("modflow6_grid" and "modflow6_model") not in dataset.attrs:
-            # TODO error
-            pass
-        else:
-            modelname = dataset.attrs["modflow6_model"].split(":")[0].lower()
-            gridtype = dataset.attrs["modflow6_grid"].lower()
-            if dis_str == "vertex":
-                if gridtype == "layered mesh":
-                    nc_dataset = DisvNetCDFMesh2d()
-            elif dis_str == "structured":
-                if gridtype == "layered mesh":
-                    nc_dataset = DisNetCDFMesh2d()
-                elif gridtype == "structured":
-                    nc_dataset = DisNetCDFStructured()
-
-        dataset.close()
-    except Exception as e:
-        print(f"Exception: {e}")
+    dataset = xr.open_dataset(fpth, engine="netcdf4")
+    if ("modflow6_grid" and "modflow6_model") not in dataset.attrs:
+        modelname = None
+        gridtype = None
+    else:
+        modelname = dataset.attrs["modflow6_model"].split(":")[0].lower()
+        gridtype = dataset.attrs["modflow6_grid"].lower()
+        if dis_str == "vertex":
+            if gridtype == "layered mesh":
+                nc_dataset = DisvNetCDFMesh2d()
+        elif dis_str == "structured":
+            if gridtype == "layered mesh":
+                nc_dataset = DisNetCDFMesh2d()
+            elif gridtype == "structured":
+                nc_dataset = DisNetCDFStructured()
+
+    dataset.close()
 
     if nc_dataset:
         fpth = Path(nc_fpth).resolve()
         nc_dataset.open(fpth)
     else:
-        raise TypeError(
+        raise Exception(
             f"Unable to load netcdf dataset for file grid "
             f'type "{gridtype}" and discretization grid '
             f'type "{dis_str}"'
@@ -380,7 +364,9 @@ def open_netcdf_dataset(nc_fpth: str, dis_type: str) -> ModelNetCDFDataset:
     return nc_dataset
 
 
-def create_netcdf_dataset(model_type, name, nc_type, nc_fname, dis):
+def create_netcdf_dataset(
+    model_type, name, nc_type, nc_fname, dis
+) -> ModelNetCDFDataset:
     assert (
         model_type.lower() == "gwf6"
         or model_type.lower() == "gwt6"
@@ -399,7 +385,7 @@ def create_netcdf_dataset(model_type, name, nc_type, nc_fname, dis):
     if nc_dataset:
         nc_dataset.create(model_type, name, nc_type, nc_fname, dis)
     else:
-        raise TypeError(
+        raise Exception(
             f"Unable to generate netcdf dataset for file type "
             f'"{nc_type.lower()}" and discretization grid type '
             f'"{dis.grid_type}"'