From 664a98b2a4f37858334912ca7e271dd98fa1fbee Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 2 Mar 2022 21:54:50 +0100
Subject: [PATCH 01/35] Add functionality for writing structured data to DISU

---
 imod/mf6/disu.py                | 267 ++++++++++++++++++++++++++++++++
 imod/mf6/pkgbase.py             |  22 ++-
 imod/templates/mf6/gwf-disu.j2  |  52 +++----
 tests/test_mf6/test_mf6_disu.py | 173 +++++++++++++++++++++
 4 files changed, 483 insertions(+), 31 deletions(-)
 create mode 100644 imod/mf6/disu.py
 create mode 100644 tests/test_mf6/test_mf6_disu.py

diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
new file mode 100644
index 000000000..da87e88e5
--- /dev/null
+++ b/imod/mf6/disu.py
@@ -0,0 +1,267 @@
+import numba as nb
+import numpy as np
+import xarray as xr
+from scipy import sparse
+
+import imod
+from imod.mf6.pkgbase import Package
+
+IntArray = np.ndarray
+
+
+@nb.njit(inline="always")
+def _number(k, i, j, nrow, ncolumn):
+    return k * (nrow * ncolumn) + i * ncolumn + j
+
+
+# @nb.njit
+def _structured_connectivity(idomain):
+    nlayer, nrow, ncolumn = idomain.shape
+    # Pre-allocate: structured connectivity implies maximum of 8 neighbors
+    nconnection = idomain.size * 8
+    ii = np.empty(nconnection, dtype=np.int32)
+    jj = np.empty(nconnection, dtype=np.int32)
+
+    connection = 0
+    for k in range(nlayer):
+        for i in range(nrow):
+            for j in range(ncolumn):
+                # Skip inactive or pass-through cells
+                if idomain[k, i, j] <= 0:
+                    continue
+
+                if j < ncolumn - 1:
+                    if idomain[k, i, j + 1] > 0:
+                        ii[connection] = _number(k, i, j, nrow, ncolumn)
+                        jj[connection] = _number(k, i, j + 1, nrow, ncolumn)
+                        connection += 1
+
+                if i < nrow - 1:
+                    if idomain[k, i + 1, j] > 0:
+                        ii[connection] = _number(k, i, j, nrow, ncolumn)
+                        jj[connection] = _number(k, i + 1, j, nrow, ncolumn)
+                        connection += 1
+
+                if k < nlayer - 1:
+                    kk = k
+                    while kk < nlayer - 1:
+                        kk += 1
+                        below = idomain[kk, i, j]
+                        if below > 0:
+                            ii[connection] = _number(k, i, j, nrow, ncolumn)
+                            jj[connection] = _number(kk, i, j, nrow, ncolumn)
+                            connection += 1
+                            break
+                        elif below == 0:
+                            break
+
+    return ii[:connection], jj[:connection]
+
+
+class UnstructuredDiscretization(Package):
+    """
+    Unstructured Discretization (DISU).
+
+    Parameters
+    ----------
+    xorigin: float
+    yorigin: float
+    top: array of floats (xr.DataArray)
+    bottom: array of floats (xr.DataArray)
+    area: array of floats (xr.DataArray)
+    iac: array of integers
+    ja: array of integers
+    ihc: array of integers
+    cl12: array of floats
+    hwva: array of floats
+    """
+
+    _pkg_id = "disu"
+    _grid_data = {
+        "top": np.float64,
+        "bottom": np.float64,
+        "area": np.float64,
+        "iac": np.int32,
+        "ja": np.int32,
+        "ihc": np.int32,
+        "cl12": np.float64,
+        "hwva": np.float64,
+    }
+    _keyword_map = {"bottom": "bot"}
+    _template = Package._initialize_template(_pkg_id)
+
+    def __init__(
+        self,
+        xorigin,
+        yorigin,
+        top,
+        bottom,
+        area,
+        iac,
+        ja,
+        ihc,
+        cl12,
+        hwva,
+    ):
+        super().__init__(locals())
+        self.dataset["xorigin"] = xorigin
+        self.dataset["yorigin"] = yorigin
+        self.dataset["top"] = top
+        self.dataset["bottom"] = bottom
+        self.dataset["area"] = area
+        self.dataset["iac"] = iac
+        self.dataset["ja"] = ja
+        self.dataset["ihc"] = ihc
+        self.dataset["cl12"] = cl12
+        self.dataset["hwva"] = hwva
+
+    def render(self, directory, pkgname, globaltimes, binary):
+        disdirectory = directory / pkgname
+        d = {}
+        d["xorigin"] = float(self.dataset["xorigin"])
+        d["yorigin"] = float(self.dataset["yorigin"])
+
+        # Dimensions
+        d["nodes"] = self.dataset["top"].size
+        d["nja"] = int(self.dataset["iac"].sum())
+
+        # Grid data
+        d["top"] = self._compose_values(
+            self.dataset["top"], disdirectory, "top", binary=binary
+        )[1][0]
+        d["bot"] = self._compose_values(
+            self["bottom"], disdirectory, "bot", binary=binary
+        )[1][0]
+        d["area"] = self._compose_values(
+            self["area"], disdirectory, "area", binary=binary
+        )[1][0]
+
+        # Connection data
+        d["iac"] = self._compose_values(
+            self["iac"], disdirectory, "iac", binary=binary
+        )[1][0]
+        d["ja"] = self._compose_values(self["ja"], disdirectory, "ja", binary=binary)[
+            1
+        ][0]
+        d["ihc"] = self._compose_values(
+            self["ihc"], disdirectory, "ihc", binary=binary
+        )[1][0]
+        d["cl12"] = self._compose_values(
+            self["cl12"], disdirectory, "cl12", binary=binary
+        )[1][0]
+        d["hwva"] = self._compose_values(
+            self["hwva"], disdirectory, "hwva", binary=binary
+        )[1][0]
+
+        return self._template.render(d)
+
+    @staticmethod
+    def from_structured(
+        top,
+        bottom,
+        idomain,
+    ):
+        x = idomain.coords["x"]
+        y = idomain.coords["y"]
+        layer = idomain.coords["layer"]
+        active = idomain.values > 0
+
+        ncolumn = x.size
+        nrow = y.size
+        nlayer = layer.size
+        size = idomain.size
+
+        dx, xmin, _ = imod.util.coord_reference(x)
+        dy, ymin, _ = imod.util.coord_reference(y)
+
+        # MODFLOW6 expects the ja values to contain the cell number first
+        # while the row should be otherwise sorted ascending.
+        # scipy.spare.csr_matrix will sort the values ascending, but
+        # would not put the cell number first. To ensure this, we use
+        # the values as well as i and j; we sort on the zeros (thereby ensuring
+        # it results as a first value per column), but the actual value
+        # is the (negative) cell number (in v).
+        ii, jj = _structured_connectivity(idomain)
+        ii += 1
+        jj += 1
+        nodes = np.arange(1, size + 1, dtype=np.int32)[active.ravel()]
+        zeros = np.zeros_like(nodes)
+        i = np.concatenate([nodes, ii, jj])
+        j = np.concatenate([zeros, jj, ii])
+        v = np.concatenate([-nodes, jj, ii])
+        csr = sparse.csr_matrix((v, (i, j)), shape=(size + 1, size + 1))
+        # The first column can be identified by its negative (node) number.
+        # This entry does not require data in ihc, cl12, hwva.
+        is_node = csr.data < 0
+
+        # Constructing the CSR matrix will have sorted all the values are
+        # required by MODFLOW6. However, we're using the original structured
+        # numbering, which includes inactive cells.
+        # For MODFLOW6, we use the reduced numbering, excluding all inactive
+        # cells. This means getting rid of empty rows (iac), generating (via
+        # cumsum) new numbers, and extracting them in the right order.
+        nnz = csr.getnnz(axis=1)
+        iac = nnz[nnz > 0]
+        ja_index = np.abs(csr.data) - 1  # Get rid of negative values temporarily.
+        ja = active.ravel().cumsum()[ja_index]
+
+        # From CSR back to COO form
+        # connectivity for every cell: n -> m
+        n = np.repeat(np.arange(size + 1), nnz) - 1
+        m = csr.indices - 1
+        # Ignore the values that do not represent n -> m connections
+        n[is_node] = 0
+        m[is_node] = 0
+
+        # Based on the row and column number differences we can derive the type
+        # of connection (unless the model is a single row or single column!).
+        diff = np.abs(n - m)
+        is_vertical = (diff > 0) & (diff % (nrow * ncolumn) == 0)  # one or more layers
+        is_x = diff == 1
+        is_y = diff == ncolumn
+        is_horizontal = is_x | is_y
+
+        # We need the indexes twice. Store for re-use.
+        # As the input is structured, we need only look at cell n, not m.
+        # (n = row, m = column off the connectivity matrix.)
+        index_x = n[is_x]
+        index_y = n[is_y]
+        index_v = n[is_vertical]
+
+        # Create flat arrays for easy indexing.
+        cellheight = top.values.ravel() - bottom.values.ravel()
+        dyy, dxx = np.meshgrid(
+            np.ones(ncolumn) * np.abs(dx),
+            np.ones(nrow) * np.abs(dy),
+            indexing="ij",
+        )
+        dyy = np.repeat(dyy, nlayer).ravel()
+        dxx = np.repeat(dxx, nlayer).ravel()
+        area = dyy * dxx
+
+        # Allocate connectivity geometry arrays, all size nja.
+        ihc = is_horizontal.astype(np.int32)
+        cl12 = np.zeros_like(ihc, dtype=np.float64)
+        hwva = np.zeros_like(ihc, dtype=np.float64)
+        # Fill.
+        cl12[is_x] = 0.5 * dxx[index_x]  # cell center to vertical face
+        cl12[is_y] = 0.5 * dyy[index_y]  # cell center to vertical face
+        cl12[is_vertical] = 0.5 * cellheight[index_v]  # cell center to horizontal face
+        hwva[is_x] = dyy[index_x]  # width
+        hwva[is_y] = dxx[index_y]  # width
+        hwva[is_vertical] = area[index_v]  # area
+
+        # Set "node" and "nja" as the dimension in accordance with MODFLOW6.
+        # Should probably be updated if we could move to UGRID3D...
+        return UnstructuredDiscretization(
+            xorigin=xmin,
+            yorigin=ymin,
+            top=xr.DataArray(top.values[active], dims=["node"]),
+            bottom=xr.DataArray(bottom.values[active], dims=["node"]),
+            area=xr.DataArray(area[active.ravel()], dims=["node"]),
+            iac=xr.DataArray(iac, dims=["node"]),
+            ja=xr.DataArray(ja, dims=["nja"]),
+            ihc=xr.DataArray(ihc, dims=["nja"]),
+            cl12=xr.DataArray(cl12, dims=["nja"]),
+            hwva=xr.DataArray(hwva, dims=["nja"]),
+        )
diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index e6dd60fa3..985cd066e 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -46,6 +46,17 @@ def disv_recarr(arrdict, layer, notnull):
     return recarr
 
 
+def disu_recarr(arrdict, layer, notnull):
+    index_spec = [("node", np.int32)]
+    field_spec = [(key, np.float64) for key in arrdict]
+    sparse_dtype = np.dtype(index_spec + field_spec)
+    # Initialize the structured array
+    nrow = notnull.sum()
+    recarr = np.empty(nrow, dtype=sparse_dtype)
+    recarr["node"] = np.argwhere(notnull) + 1
+    return recarr
+
+
 class Package(abc.ABC):
     """
     Package is used to share methods for specific packages with no time
@@ -206,7 +217,11 @@ def to_sparse(self, arrdict, layer):
         notnull = ~np.isnan(data)
 
         if isinstance(self.dataset, xr.Dataset):
-            recarr = dis_recarr(arrdict, layer, notnull)
+            # TODO
+            if ("node" in self.dataset.dims) or ("nja" in self.dataset.dims):
+                recarr = disu_recarr(arrdict, layer, notnull)
+            else:
+                recarr = dis_recarr(arrdict, layer, notnull)
         elif isinstance(self.dataset, xu.UgridDataset):
             recarr = disv_recarr(arrdict, layer, notnull)
         else:
@@ -287,7 +302,10 @@ def render(self, directory, pkgname, globaltimes, binary):
     @staticmethod
     def _is_xy_data(obj):
         if isinstance(obj, (xr.DataArray, xr.Dataset)):
-            xy = "x" in obj.dims and "y" in obj.dims
+            # TODO
+            xy = ("x" in obj.dims and "y" in obj.dims) or (
+                "node" in obj.dims or "nja" in obj.dims
+            )
         elif isinstance(obj, (xu.UgridDataArray, xu.UgridDataset)):
             xy = obj.ugrid.grid.face_dimension in obj.dims
         else:
diff --git a/imod/templates/mf6/gwf-disu.j2 b/imod/templates/mf6/gwf-disu.j2
index c00936b55..f1c2d7576 100644
--- a/imod/templates/mf6/gwf-disu.j2
+++ b/imod/templates/mf6/gwf-disu.j2
@@ -1,49 +1,43 @@
 begin options
-{%- if length_units is defined -%}  length_units {{length_units}}{%- endif -%}
-{%- if nogrb is defined -%}  nogrb{%- endif -%}
-{%- if xorigin is defined -%}  xorigin {{xorigin}}{%- endif -%}
-{%- if yorigin is defined -%}  yorigin {{yorigin}}{%- endif -%}
-{%- if angrot is defined -%}  angrot {{angrot}}{%- endif -%}
+{% if length_units is defined %}  length_units {{length_units}}
+{% endif %}
+{%- if nogrb is defined %}  nogrb
+{% endif %}
+{%- if xorigin is defined %}  xorigin {{xorigin}}
+{% endif %}
+{%- if yorigin is defined %}  yorigin {{yorigin}}
+{% endif %}
+{%- if angrot is defined %}  angrot {{angrot}}
+{% endif -%}
 end options
 
 begin dimensions
   nodes {{nodes}}
   nja {{nja}}
-{%- if nvert is defined -%}  nvert {{nvert}}{%- endif -%}
+{%- if nvert is defined %}  nvert {{nvert}}
+{% endif %}
 end dimensions
 
 begin griddata
   top
-    {top}
+    {{top}}
   bot
-    {bot}
+    {{bot}}
   area
-    {area}
+    {{area}}
 end griddata
 
 begin connectiondata
   iac
-    {iac}
+    {{iac}}
   ja
-    {ja}
+    {{ja}}
   ihc
-    {ihc}
+    {{ihc}}
   cl12
-    {cl12}
+    {{cl12}}
   hwva
-    {hwva}
-{%- if angldegx is defined -%}  angldegx
-    {angldegx}{%- endif -%}
-end connectiondata
-
-begin vertices
-  {{iv}} {{xv}} {{yv}}
-  {{iv}} {{xv}} {{yv}}
-  ...
-end vertices
-
-begin cell2d
-  {{icell2d}} {{xc}} {{yc}} {{ncvert}} {{icvert(ncvert)}}
-  {{icell2d}} {{xc}} {{yc}} {{ncvert}} {{icvert(ncvert)}}
-  ...
-end cell2d
+    {{hwva}}
+{%- if angldegx is defined %}  angldegx
+    {angldegx}{% endif %}
+end connectiondata
\ No newline at end of file
diff --git a/tests/test_mf6/test_mf6_disu.py b/tests/test_mf6/test_mf6_disu.py
new file mode 100644
index 000000000..ca552f9b1
--- /dev/null
+++ b/tests/test_mf6/test_mf6_disu.py
@@ -0,0 +1,173 @@
+import pathlib
+import textwrap
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from imod.mf6 import disu
+
+
+@pytest.fixture(scope="function")
+def structured_dis():
+    coords = {
+        "layer": [1, 2, 3],
+        "y": [25.0, 15.0, 5.0],
+        "x": [50.0, 70.0, 90.0, 110.0],
+    }
+    dims = ["layer", "y", "x"]
+    idomain = xr.DataArray(np.ones((3, 3, 4), dtype=np.int32), coords, dims)
+    top = xr.DataArray(np.ones((3, 3, 4)), coords, dims)
+    bottom = xr.DataArray(np.ones((3, 3, 4)), coords, dims)
+    top.values[0] = 45.0
+    top.values[1] = 30.0
+    top.values[2] = 15.0
+    bottom.values[0] = 30.0
+    bottom.values[1] = 15.0
+    bottom.values[2] = 0.0
+    return idomain, top, bottom
+
+
+def connectivity_checks(dis):
+    def check_symmetry(i, j, v):
+        order0 = np.lexsort((j, i))
+        order1 = np.lexsort((i, j))
+        assert np.allclose(v[order0], v[order1])
+
+    ncell = dis.dataset["node"].size
+    i = np.repeat(np.arange(1, ncell + 1), dis.dataset["iac"].values) - 1
+    j = dis.dataset["ja"].values - 1
+    diff = np.abs(i - j)
+    ihc = dis.dataset["ihc"].values
+    hwva = dis.dataset["hwva"].values
+    cl12 = dis.dataset["cl12"].values
+
+    node = i == j
+    connection = i != j
+    vertical = (ihc == 0) & connection
+
+    assert dis.dataset["iac"].values.sum() == j.size
+    assert i.min() == 0
+    assert j.min() == 0
+    assert i.max() == ncell - 1
+    assert j.max() == ncell - 1
+    assert (diff[node] == 0).all()
+    assert (diff[connection] > 0).all()
+    assert (cl12[node] == 0).all()
+    assert (cl12[connection] != 0).all()
+    assert (hwva[node] == 0).all()
+    assert (hwva[connection] != 0).all()
+    assert (diff[vertical] > 1).all()
+    assert np.allclose(cl12[vertical], 7.5)
+    assert np.allclose(hwva[vertical], 200.0)
+    check_symmetry(i, j, hwva)
+
+
+def test_cell_number():
+    nrow = 4
+    ncol = 3
+    assert disu._number(0, 0, 0, nrow, ncol) == 0
+    assert disu._number(0, 0, 1, nrow, ncol) == 1
+    assert disu._number(0, 1, 0, nrow, ncol) == 3
+    assert disu._number(0, 1, 1, nrow, ncol) == 4
+    assert disu._number(1, 0, 0, nrow, ncol) == 12
+    assert disu._number(1, 0, 1, nrow, ncol) == 13
+    assert disu._number(1, 1, 0, nrow, ncol) == 15
+    assert disu._number(1, 1, 1, nrow, ncol) == 16
+
+
+def test_structured_connectivity_full():
+    idomain = np.ones((3, 3, 4), dtype=np.int32)
+    i, j = disu._structured_connectivity(idomain)
+
+    expected_i = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3])
+    expected_j = np.array([1, 4, 12, 2, 5, 13, 3, 6, 14, 7, 15])
+    assert np.array_equal(i[:11], expected_i)
+    assert np.array_equal(j[:11], expected_j)
+
+    idomain[1, 0, 1] = -1
+    idomain[1, 0, 2] = -1
+    i, j = disu._structured_connectivity(idomain)
+
+    expected_i = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3])
+    expected_j = np.array([1, 4, 12, 2, 5, 25, 3, 6, 26, 7, 15])
+    assert np.array_equal(i[:11], expected_i)
+    assert np.array_equal(j[:11], expected_j)
+
+
+def test_from_structured(structured_dis):
+    idomain, top, bottom = structured_dis
+    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    assert np.allclose(dis.dataset["xorigin"], 40.0)
+    assert np.allclose(dis.dataset["yorigin"], 0.0)
+    connectivity_checks(dis)
+
+    # Now disable some cells, create one pass-through
+    idomain.values[1, 0, 1] = -1
+    idomain.values[:, 0, 0] = 0
+    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    connectivity_checks(dis)
+
+
+def test_render(structured_dis):
+    idomain, top, bottom = structured_dis
+    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+
+    directory = pathlib.Path("mymodel")
+    actual = dis.render(directory, "disu", None, True)
+
+    expected = textwrap.dedent(
+        """\
+        begin options
+          xorigin 40.0
+          yorigin 0.0
+        end options
+
+        begin dimensions
+          nodes 36
+          nja 186
+        end dimensions
+
+        begin griddata
+          top
+            open/close mymodel/disu/top.bin (binary)
+          bot
+            open/close mymodel/disu/bot.bin (binary)
+          area
+            open/close mymodel/disu/area.bin (binary)
+        end griddata
+
+        begin connectiondata
+          iac
+            open/close mymodel/disu/iac.bin (binary)
+          ja
+            open/close mymodel/disu/ja.bin (binary)
+          ihc
+            open/close mymodel/disu/ihc.bin (binary)
+          cl12
+            open/close mymodel/disu/cl12.bin (binary)
+          hwva
+            open/close mymodel/disu/hwva.bin (binary)
+        end connectiondata"""
+    )
+    assert actual == expected
+
+
+def test_write(structured_dis, tmp_path):
+    idomain, top, bottom = structured_dis
+    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis.write(tmp_path, "disu", None, True)
+
+    assert (tmp_path / "disu.disu").exists
+    names = [
+        "top.bin",
+        "bot.bin",
+        "area.bin",
+        "iac.bin",
+        "ja.bin",
+        "ihc.bin",
+        "cl12.bin",
+        "hwva.bin",
+    ]
+    for name in names:
+        assert (tmp_path / name).exists

From 97135b2ab9653f4536e091e6e59b35f6f4baed6d Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 3 Mar 2022 15:43:18 +0100
Subject: [PATCH 02/35] Rename to LowLevel; add pkg conversion to disu

---
 imod/mf6/disu.py                |  4 ++--
 imod/mf6/pkgbase.py             | 20 ++++++++++++++++++--
 tests/test_mf6/test_mf6_disu.py |  8 ++++----
 3 files changed, 24 insertions(+), 8 deletions(-)

diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
index da87e88e5..68db74130 100644
--- a/imod/mf6/disu.py
+++ b/imod/mf6/disu.py
@@ -58,7 +58,7 @@ def _structured_connectivity(idomain):
     return ii[:connection], jj[:connection]
 
 
-class UnstructuredDiscretization(Package):
+class LowLevelUnstructuredDiscretization(Package):
     """
     Unstructured Discretization (DISU).
 
@@ -253,7 +253,7 @@ def from_structured(
 
         # Set "node" and "nja" as the dimension in accordance with MODFLOW6.
         # Should probably be updated if we could move to UGRID3D...
-        return UnstructuredDiscretization(
+        return LowLevelUnstructuredDiscretization(
             xorigin=xmin,
             yorigin=ymin,
             top=xr.DataArray(top.values[active], dims=["node"]),
diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index 985cd066e..f5b7afbe5 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -416,7 +416,18 @@ def write_netcdf(self, directory, pkgname, aggregate_layers=False):
 
         spatial_ds.to_netcdf(path)
         return has_dims
-
+        
+    def _to_disu(self, numbers):
+        structured = self.dataset.expand_dims("layer")
+        # Stack will automatically broadcast to layer
+        dataset = structured.stack(node=("layer", "y", "x"))
+        layers = structured["layer"].values
+        ncell_per_layer = structured["y"].size * structured["x"].size
+        offset = (layers - 1) * ncell_per_layer
+        index = np.add.outer(offset, np.arange(ncell_per_layer))
+        dataset = dataset.assign_coords(node=numbers[index])
+        return self.__class__(**dataset)
+ 
 
 class BoundaryCondition(Package, abc.ABC):
     """
@@ -455,10 +466,15 @@ def write_datafile(self, outpath, ds, binary):
         """
         Writes a modflow6 binary data file
         """
-        layer = ds["layer"].values
+        if "layer" in ds:
+            layer = ds["layer"].values
+        else:
+            layer = None
+
         arrdict = self._ds_to_arrdict(ds)
         sparse_data = self.to_sparse(arrdict, layer)
         outpath.parent.mkdir(exist_ok=True, parents=True)
+
         if binary:
             self._write_binaryfile(outpath, sparse_data)
         else:
diff --git a/tests/test_mf6/test_mf6_disu.py b/tests/test_mf6/test_mf6_disu.py
index ca552f9b1..ed6075540 100644
--- a/tests/test_mf6/test_mf6_disu.py
+++ b/tests/test_mf6/test_mf6_disu.py
@@ -97,7 +97,7 @@ def test_structured_connectivity_full():
 
 def test_from_structured(structured_dis):
     idomain, top, bottom = structured_dis
-    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
     assert np.allclose(dis.dataset["xorigin"], 40.0)
     assert np.allclose(dis.dataset["yorigin"], 0.0)
     connectivity_checks(dis)
@@ -105,13 +105,13 @@ def test_from_structured(structured_dis):
     # Now disable some cells, create one pass-through
     idomain.values[1, 0, 1] = -1
     idomain.values[:, 0, 0] = 0
-    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
     connectivity_checks(dis)
 
 
 def test_render(structured_dis):
     idomain, top, bottom = structured_dis
-    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
 
     directory = pathlib.Path("mymodel")
     actual = dis.render(directory, "disu", None, True)
@@ -155,7 +155,7 @@ def test_render(structured_dis):
 
 def test_write(structured_dis, tmp_path):
     idomain, top, bottom = structured_dis
-    dis = disu.UnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
     dis.write(tmp_path, "disu", None, True)
 
     assert (tmp_path / "disu.disu").exists

From fd5e74359ee890991190490b4ba3ee8c1a7f416d Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 4 Apr 2022 15:56:37 +0200
Subject: [PATCH 03/35] Make a start with reading MODFLOW6 input files

---
 imod/mf6/dis.py                   |  34 ++++
 imod/mf6/model.py                 |  99 ++++++++++
 imod/mf6/pkgbase.py               |  29 +++
 imod/mf6/read_input/__init__.py   | 298 ++++++++++++++++++++++++++++++
 imod/mf6/read_input/common.py     | 152 +++++++++++++++
 imod/mf6/read_input/grid_data.py  | 197 ++++++++++++++++++++
 imod/mf6/read_input/list_input.py | 102 ++++++++++
 imod/mf6/simulation.py            |  36 ++++
 8 files changed, 947 insertions(+)
 create mode 100644 imod/mf6/read_input/__init__.py
 create mode 100644 imod/mf6/read_input/common.py
 create mode 100644 imod/mf6/read_input/grid_data.py
 create mode 100644 imod/mf6/read_input/list_input.py

diff --git a/imod/mf6/dis.py b/imod/mf6/dis.py
index 83250f834..0998a9427 100644
--- a/imod/mf6/dis.py
+++ b/imod/mf6/dis.py
@@ -1,8 +1,13 @@
+from pathlib import Path
+
 import numpy as np
+import xarray as xr
 
 import imod
 from imod.mf6.pkgbase import Package
 
+from .read_input import read_dis_blockfile
+
 
 class StructuredDiscretization(Package):
     """
@@ -53,6 +58,35 @@ def _delrc(self, dx):
         else:
             raise ValueError(f"Unhandled type of {dx}")
 
+    @classmethod
+    def open(cls, path: str, simroot: Path) -> "StructuredDiscretization":
+        content = read_dis_blockfile(simroot / path, simroot)
+        nlay = content["nlay"]
+
+        griddata = content["griddata"]
+        xmin = float(content.get("xorigin", 0.0))
+        ymin = float(content.get("yorigin", 0.0))
+        dx = griddata.pop("delr")
+        dy = griddata.pop("delc")
+        coords = {
+            "layer": np.arange(1, nlay + 1),
+            "y": ((ymin + np.cumsum(dy)) - 0.5 * dx)[::-1],
+            "x": (xmin + np.cumsum(dx)) - 0.5 * dx,
+        }
+        dims = ("layer", "y", "x")
+
+        inverse_keywords = {v: k for k, v in cls._keyword_map.items()}
+        variables = {}
+        for key, value in griddata.items():
+            invkey = inverse_keywords.get(key, key)
+            variables[invkey] = xr.DataArray(
+                value,
+                coords,
+                dims,
+            )
+
+        return cls(**variables)
+
     def render(self, directory, pkgname, globaltimes, binary):
         disdirectory = directory / pkgname
         d = {}
diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 2db089847..4208fb61d 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -1,12 +1,16 @@
 import collections
 import pathlib
+from typing import Dict, Type, Union
 
 import cftime
 import jinja2
 import numpy as np
 
+from imod import mf6
 from imod.mf6 import qgs_util
 
+from .read_input import read_gwf_namefile
+
 
 class Model(collections.UserDict):
     def __setitem__(self, key, value):
@@ -25,6 +29,31 @@ class GroundwaterFlowModel(Model):
 
     _pkg_id = "model"
 
+    @staticmethod
+    def _PACKAGE_CLASSES() -> Dict[str, Type]:
+        return {
+            package._pkg_id: package
+            for package in (
+                mf6.ConstantHead,
+                mf6.StructuredDiscretization,
+                mf6.VerticesDiscretization,
+                mf6.Drainage,
+                mf6.Evapotranspiration,
+                mf6.GeneralHeadBoundary,
+                mf6.InitialConditions,
+                mf6.Solution,
+                mf6.NodePropertyFlow,
+                mf6.OutputControl,
+                mf6.Recharge,
+                mf6.River,
+                mf6.SpecificStorage,
+                mf6.Storage,
+                mf6.StorageCoefficient,
+                mf6.WellDisStructured,
+                mf6.WellDisVertices,
+            )
+        }
+
     def _initialize_template(self):
         loader = jinja2.PackageLoader("imod", "templates/mf6")
         env = jinja2.Environment(loader=loader, keep_trailing_newline=True)
@@ -109,6 +138,76 @@ def render(self, modelname):
         d["packages"] = packages
         return self._template.render(d)
 
+    @classmethod
+    def open(
+        cls,
+        path: Union[pathlib.Path, str],
+        simroot: pathlib.Path,
+        globaltimes: np.ndarray,
+    ):
+        content = read_gwf_namefile(simroot / path)
+        model = cls(**content)
+
+        # Search for the DIS/DISV/DISU package first. This provides us with
+        # the coordinates and dimensions to instantiate the other packages.
+        classes = cls._PACKAGE_CLASSES()
+        dis_packages = [
+            tup for tup in content["packages"] if tup[0] in ("dis6", "disv6", "disu6")
+        ]
+        packages = [
+            tup
+            for tup in content["packages"]
+            if tup[0] not in ("dis6", "disv6", "disu6")
+        ]
+
+        if len(dis_packages) > 1:
+            raise ValueError(f"More than one DIS/DISV/DISU package in {path}")
+
+        disftype, disfname, dispname = dis_packages[0]
+        diskey = disftype[:-1]
+        package = classes[diskey]
+        if dispname is None:
+            disname = diskey
+
+        dis_package = package.open(
+            disfname,
+            simroot,
+        )
+        model[disname] = dis_package
+        shape = dis_package["idomain"].shape
+        coords = dis_package["idomain"].coords
+        dims = dis_package["idomain"].dims
+
+        # Now read the rest of the packages.
+        seen = set()
+        for i, (ftype, fname, pname) in packages:
+
+            key = ftype[:-1]  # Remove the last number (riv6 -> riv).
+            package = classes[key]
+
+            # Fill in a name if none is given.
+            if pname is None:
+                pname = key
+
+            # Ensure a unique name is generated.
+            if pname in seen:
+                pkgname = f"{pname}_{i+1}"
+            else:
+                pkgname = pname
+            seen.add(pkgname)
+
+            # Create the package and add it to the model.
+            model[pkgname] = package.open(
+                fname,
+                simroot,
+                shape,
+                coords,
+                dims,
+                globaltimes,
+            )
+
+        return model
+
     def write(self, directory, modelname, globaltimes, binary=True):
         """
         Write model namefile
diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index e6dd60fa3..ef43b11af 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -6,6 +6,8 @@
 import xarray as xr
 import xugrid as xu
 
+from .read_input import get_function_kwargnames, read_blockfile
+
 
 def dis_recarr(arrdict, layer, notnull):
     # Define the numpy structured array dtype
@@ -525,6 +527,33 @@ def write(self, directory, pkgname, globaltimes, binary):
             binary=binary,
         )
 
+    @classmethod
+    def open(cls, path, simroot, shape, coords, dims, globaltimes):
+        content = read_blockfile(
+            simroot / path,
+            simroot,
+            fields=cls._period_data,
+            shape=shape,
+            dims=dims,
+        )
+
+        period_index = content.pop("period_index")
+        period_data = content.pop("period_data")
+        coords = coords.copy()
+        coords["time"] = globaltimes[period_index]
+        dims = ("time",) + dims
+
+        for field, data in period_data.items():
+            content[field] = xr.DataArray(data, coords, dims)
+
+        filtered_content = {
+            k: v
+            for k, v in content.items()
+            if k in get_function_kwargnames(cls.__init__)
+        }
+
+        return cls(**filtered_content)
+
 
 class AdvancedBoundaryCondition(BoundaryCondition, abc.ABC):
     """
diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
new file mode 100644
index 000000000..cb6e3e3b9
--- /dev/null
+++ b/imod/mf6/read_input/__init__.py
@@ -0,0 +1,298 @@
+"""
+Utilities for reading MODFLOW6 input files.
+"""
+import inspect
+import warnings
+from collections import defaultdict
+from pathlib import Path
+from typing import IO, Any, Callable, Dict, List, Tuple
+
+import dask.array
+import numpy as np
+
+from .common import (
+    advance_to_header,
+    advance_to_period,
+    parse_dimension,
+    parse_option,
+    read_iterable_block,
+    read_key_value_block,
+)
+from .grid_data import read_griddata, shape_to_max_rows
+from .list_input import read_listinput
+
+
+def parse_model(stripped: str, fname: str) -> Tuple[str, str, str]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 3:
+        return separated
+    else:
+        raise ValueError(
+            "ftype, fname and pname expected. Received instead: "
+            f"{','.join(separated)} in file {fname}"
+        )
+
+
+def parse_exchange(stripped: str, fname: str) -> Tuple[str, str, str, str]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 4:
+        return separated
+    else:
+        raise ValueError(
+            "exgtype, exgfile, exgmnamea, exgmnameb expected. Received instead: "
+            f"{','.join(separated)} in file {fname}"
+        )
+
+
+def parse_solutiongroup(stripped: str, fname: str) -> Tuple[str, str]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords > 2:
+        return separated[0], separated[1], separated[2:]
+    else:
+        raise ValueError(
+            "Expected at least three entries: slntype, slnfname, and one model name. "
+            f"Received instead: {','.join(separated)} in file {fname}"
+        )
+
+
+def parse_package(stripped: str, fname: str) -> Tuple[str, str, str]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 2:
+        ftype, fname = separated
+        pname = ftype
+    elif nwords == 3:
+        ftype, fname, pname = separated
+    else:
+        raise ValueError(
+            "Expected ftype, fname, and optionally pname. "
+            f"Received instead: {','.join(separated)} in file {fname}"
+        )
+    return ftype, fname, pname
+
+
+def parse_tdis_perioddata(stripped: str, fname: str) -> Tuple[float, int, float]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 3:
+        return float(separated[0]), int(separated[1]), float(separated[2])
+    else:
+        raise ValueError(
+            "perlen, nstp, tsmult expected. Received instead: "
+            f"{','.join(separated)} in file {fname}"
+        )
+
+
+def read_tdis(path: Path) -> Dict[str, Any]:
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "dimensions")
+        dimensions = read_key_value_block(f, parse_dimension)
+        advance_to_header(f, "perioddata")
+        content["perioddata"] = read_iterable_block(f, parse_tdis_perioddata)
+    return {**content, **dimensions}
+
+
+def read_dis_blockfile(path: Path, simroot: Path) -> Dict[str, Any]:
+    def delr_max_rows(shape, layered) -> Tuple[int, int]:
+        if layered:
+            raise ValueError(f"DELR section in {path} is LAYERED")
+        return 1, shape[2]
+
+    def delc_max_rows(shape, layered) -> int:
+        if layered:
+            raise ValueError(f"DELC section in {path} is LAYERED")
+        return 1, shape[1]
+
+    sections = {
+        "delr": (np.float64, delr_max_rows),
+        "delc": (np.float64, delc_max_rows),
+        "top": (np.float64, shape_to_max_rows),
+        "botm": (np.float64, shape_to_max_rows),
+        "idomain": (np.int32, shape_to_max_rows),
+    }
+
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "dimensions")
+        dimensions = read_key_value_block(f, parse_dimension)
+        shape = (dimensions["nlay"], dimensions["nrow"], dimensions["ncol"])
+        advance_to_header(f, "griddata")
+        content["griddata"] = read_griddata(f, simroot, sections, shape)
+
+    return {**content, **dimensions}
+
+
+def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
+    # https://numpy.org/doc/stable/reference/arrays.datetime.html#datetime-units
+    TIME_UNITS = {
+        "unknown": None,
+        "seconds": "s",
+        "minutes": "m",
+        "hours": "h",
+        "days": "D",
+        "years": "Y",
+    }
+    unit = TIME_UNITS.get(tdis["units"])
+
+    start = None
+    if "start_date_time" in tdis:
+        try:
+            start = np.datetime64(tdis["start_date_time"])
+        except ValueError:
+            pass
+
+    cumulative_length = np.cumsum([entry[0] for entry in tdis["perioddata"]])
+    if unit is not None and start is not None:
+        timedelta = np.timedelta64(cumulative_length, unit)
+        times = np.full(timedelta.size, start)
+        times[1:] += timedelta[:-1]
+    else:
+        possibilities = ",".join(list(TIME_UNITS.keys())[1:])
+        warnings.warn(
+            "Cannot convert time coordinate to datetime. "
+            "Falling back to (unitless) floating point time coordinate. "
+            f"time_units should be one of: {possibilities}; "
+            "start_date_time should be according to ISO 8601."
+        )
+        times = np.full(timedelta.size, 0.0)
+        times[1:] += cumulative_length[:-1]
+
+    return times
+
+
+def read_solver(path: Path) -> Dict[str, str]:
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        options = read_key_value_block(f, parse_option)
+        advance_to_header(f, "nonlinear")
+        nonlinear = read_key_value_block(f, parse_option)
+        advance_to_header(f, "linear")
+        linear = read_key_value_block(f, parse_option)
+    return {**options, **nonlinear, **linear}
+
+
+def read_sim_namefile(path: Path) -> Dict[str, str]:
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "timing")
+        timing = read_key_value_block(f, parse_option)
+        advance_to_header(f, "models")
+        content["models"] = read_iterable_block(f, parse_model)
+        advance_to_header(f, "exchanges")
+        content["exchanges"] = read_iterable_block(f, parse_exchange)
+        advance_to_header(f, "solutiongroup 1")
+        content["solutiongroup 1"] = read_iterable_block(f, parse_solutiongroup)
+    return {**content, **timing}
+
+
+def read_gwf_namefile(path: Path) -> Dict[str, Any]:
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "packages")
+        content["packages"] = read_iterable_block(f, parse_package)
+    return content
+
+
+def read_package_periods(
+    f: IO[str],
+    simroot: Path,
+    dtype: type,
+    fields: Tuple[str],
+    index_columns: Tuple[str],
+    max_rows: int,
+    shape: Tuple[int],
+) -> Tuple[List[int], Dict[str, dask.array.Array]]:
+    """
+    Read blockfile periods section of a "standard" MODFLOW6 package: RIV, DRN,
+    GHB, etc.
+
+    Parameters
+    ----------
+    f: IO[str]
+        File handle.
+    simroot:
+    """
+    periods = []
+    dask_lists = defaultdict(list)
+    key = advance_to_period(f)
+
+    while key is not None:
+        # Read the recarrays, already separated into dense arrays.
+        variable_values = read_listinput(
+            f=f,
+            simroot=simroot,
+            dtype=dtype,
+            fields=fields,
+            index_columns=index_columns,
+            max_rows=max_rows,
+            shape=shape,
+        )
+
+        # Group them by field (e.g. cond, head, etc.)
+        for field, values in zip(fields, variable_values):
+            dask_lists[field].append(values)
+
+        # Store number and advance to next period
+        periods.append(key)
+        key = advance_to_period(f)
+
+    # Create a dictionary of arrays
+    variable_values = {
+        field: dask.array.stack(dask_list, axis=0)
+        for field, dask_list in dask_lists.items()
+    }
+    return periods, variable_values
+
+
+def read_blockfile(
+    path: Path,
+    simroot: Path,
+    fields: Tuple[str],
+    shape: Tuple[int],
+    dims: Tuple[str],
+) -> Dict[str, Any]:
+    """
+    Read blockfile of a "standard" MODFLOW6 package: RIV, DRN, GHB, etc.
+    """
+    ndim = len(dims)
+    index_columns = {
+        1: ("node",),
+        2: ("layer", "cell2d"),
+        3: ("layer", "row", "column"),
+    }.get(ndim)
+    if index_columns is None:
+        raise ValueError(f"Length of dimensions should be 1, 2, or 3. Received: {ndim}")
+
+    dtype = np.dtype(
+        [(column, np.int32) for column in index_columns]
+        + [(field, np.float64) for field in fields]
+    )
+
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "dimensions")
+        dimensions = read_key_value_block(f, parse_dimension)
+        content["period_index"], content["period_data"] = read_package_periods(
+            f=f,
+            simroot=simroot,
+            dtype=dtype,
+            fields=fields,
+            index_columns=index_columns,
+            max_rows=dimensions["maxbound"],
+            shape=shape,
+        )
+
+    return content
+
+
+def get_function_kwargnames(f: Callable) -> List[str]:
+    return inspect.getfullargspec(f).args
diff --git a/imod/mf6/read_input/common.py b/imod/mf6/read_input/common.py
new file mode 100644
index 000000000..24d8b1cfb
--- /dev/null
+++ b/imod/mf6/read_input/common.py
@@ -0,0 +1,152 @@
+"""
+Commonly shared utilities for reading MODFLOW6 input files.
+"""
+from pathlib import Path
+from typing import IO, Any, Callable, Dict, List, Tuple
+
+import numpy as np
+
+
+def end_of_file(line: str) -> bool:
+    return line == ""
+
+
+def strip_line(line: str) -> str:
+    # remove possible comments
+    before, _, _ = line.partition("#")
+    return before.strip().lower()
+
+
+def find_entry(line: str, entry: str, type: type) -> str:
+    if entry not in line:
+        return None
+    else:
+        _, after = line.split(entry)
+        return type(after.split()[0])
+
+
+def read_internal(f: IO[str], max_rows: int, dtype: type) -> np.ndarray:
+    return np.loadtxt(
+        fname=f,
+        dtype=dtype,
+        max_rows=max_rows,
+    )
+
+
+def read_external_binaryfile(path: Path, dtype: type) -> np.ndarray:
+    return np.fromfile(
+        file=path,
+        dtype=dtype,
+        sep="",
+    )
+
+
+def read_external_textfile(path: Path, dtype: type) -> np.ndarray:
+    return np.loadtxt(
+        fname=path,
+        dtype=dtype,
+    )
+
+
+def advance_to_header(f: IO[str], section) -> None:
+    line = None
+    start = f"begin {section}"
+    # Empty line is end-of-file
+    while not end_of_file(line):
+        line = f.readline()
+        stripped = line.strip().lower()
+        # Return if start has been found
+        if stripped == start:
+            return
+        # Continue if line is empty
+        elif stripped == "":
+            continue
+        # Raise if the line has (unexpected) content
+        else:
+            break
+    # Also raise if no further content is found
+    raise ValueError(f'"{start}" is not present in file {f.name}')
+
+
+def read_key_value_block(f: IO[str], parse: Callable) -> Dict[str, str]:
+    fname = f.name
+    content = {}
+    line = None
+    while not end_of_file(line):
+        line = f.readline()
+        stripped = strip_line(line)
+        # Return if end has been found
+        if stripped[:3] == "end":
+            return content
+        # Continue in case of an empty line
+        elif stripped == "":
+            continue
+        # Valid entry
+        else:
+            key, value = parse(stripped, fname)
+            content[key] = value
+
+    # Also raise if no further content is found
+    raise ValueError(f'"end" of block is not present in file {fname}')
+
+
+def read_iterable_block(f: IO[str], parse: Callable) -> List[Any]:
+    fname = f.name
+    content = []
+    line = None
+    while not end_of_file(line):
+        line = f.readline()
+        stripped = strip_line(line)
+        # Return if end has been found
+        if stripped[:3] == "end":
+            return content
+        # Continue in case of an empty line
+        elif stripped == "":
+            continue
+        # Valid entry
+        else:
+            content.append(parse(stripped, fname))
+
+    # Also raise if no further content is found
+    raise ValueError(f'"end" of block is not present in file {fname}')
+
+
+def parse_option(stripped: str, fname: str) -> Tuple[str, str]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 1:
+        key = separated[0]
+        value = True
+    elif nwords == 2:
+        key, value = separated
+    else:
+        raise ValueError(
+            "More than two words found in block:"
+            f"{','.join(separated)} in file {fname}"
+        )
+    return key, value
+
+
+def parse_dimension(stripped: str, fname: str) -> Tuple[str, int]:
+    key, value = parse_option(stripped, fname)
+    return key, int(value)
+
+
+def advance_to_period(f: IO[str]) -> int:
+    line = None
+    start = "begin period"
+    # Empty line is end-of-file
+    while not end_of_file(line):
+        line = f.readline()
+        stripped = strip_line(line)
+        # Return if start has been found
+        if stripped[:12] == start:
+            return int(stripped.split()[2])
+        # Continue if line is empty
+        elif stripped == "":
+            continue
+        # Raise if the line has (unexpected) content
+        else:
+            break
+    else:
+        return None
diff --git a/imod/mf6/read_input/grid_data.py b/imod/mf6/read_input/grid_data.py
new file mode 100644
index 000000000..9623d62d6
--- /dev/null
+++ b/imod/mf6/read_input/grid_data.py
@@ -0,0 +1,197 @@
+"""
+Utilities for reading MODFLOW6 GRID DATA input.
+"""
+from pathlib import Path
+from typing import IO, Any, Callable, Dict, Tuple
+
+import dask.array
+import numpy as np
+
+from .common import (
+    end_of_file,
+    find_entry,
+    read_external_binaryfile,
+    read_external_textfile,
+    read_internal,
+    strip_line,
+)
+
+
+def advance_to_griddata_section(f: IO[str], section: str) -> None:
+    line = None
+    # Empty line is end-of-file
+    while not end_of_file(line):
+        line = f.readline()
+        stripped = strip_line(line)
+        # Return if start has been found
+        if stripped == section:
+            return False
+        elif stripped == f"{section} layered":
+            return True
+        # Continue if line is empty
+        elif stripped == "":
+            continue
+        # Raise if the line has (unexpected) content
+        else:
+            break
+    raise ValueError(f'"{section}" is not present in file {f.name}')
+
+
+def shape_to_max_rows(shape: Tuple[int], layered: bool) -> int:
+    ndim = len(shape)
+    if ndim == 3:
+        nlayer, nrow, _ = shape
+        max_rows_layered = nrow
+        max_rows = nlayer * nrow
+
+    elif ndim == 2:
+        nlayer, _ = shape
+        max_rows_layered = 1
+        max_rows = nlayer
+
+    elif ndim == 1:
+        if layered:
+            raise ValueError(
+                "LAYERED section detected. DISU does not support LAYERED input."
+            )
+        nlayer = None
+        max_rows_layered = None
+        max_rows = 1
+
+    else:
+        raise ValueError(
+            f"length of shape should be 1, 2, or 3. Received shape of length: {ndim}"
+        )
+
+    if layered:
+        return max_rows_layered, shape[1:]
+    else:
+        return max_rows, shape
+
+
+def constant(value: Any, shape: Tuple[int], dtype: type) -> dask.array.Array:
+    return dask.array.full(shape, value, dtype=dtype)
+
+
+def read_internal_griddata(
+    f: IO[str], dtype: type, shape: Tuple[int], max_rows: int
+) -> np.ndarray:
+    return read_internal(f, max_rows, dtype).reshape(shape)
+
+
+def read_external_griddata(
+    path: Path, dtype: type, shape: Tuple[int], binary: bool
+) -> np.ndarray:
+    if binary:
+        a = read_external_binaryfile(path, dtype)
+    else:
+        a = read_external_textfile(path, dtype)
+    return a.reshape(shape)
+
+
+def read_array(
+    f: IO[str],
+    simroot: Path,
+    dtype: type,
+    max_rows: int,
+    shape: Tuple[int],
+) -> dask.array.Array:
+    """
+    MODFLOW6 READARRAY functionality for grid data.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    array: dask.array.Array
+    """
+    firstline = f.readline()
+    stripped = strip_line(firstline)
+    separated = stripped.split()
+    first = separated[0]
+
+    if first == "constant":
+        factor = None
+        array = constant(separated[1], shape, dtype)
+    elif first == "internal":
+        factor = find_entry(stripped, "factor", float)
+        a = read_internal_griddata(f, dtype, shape, max_rows)
+        array = dask.array.from_array(a)
+    elif first == "open/close":
+        factor = find_entry(stripped, "factor", float)
+        fname = separated[1]
+        binary = "(binary)" in stripped
+        path = simroot / fname
+        a = dask.delayed(read_external_griddata)(path, dtype, shape, binary)
+        array = dask.array.from_delayed(a, shape=shape, dtype=dtype)
+    else:
+        raise ValueError(
+            'Expected "constant", "internal" or "open/close". Received instead: '
+            f"{stripped}"
+        )
+
+    if factor is not None:
+        array = array * factor
+
+    return array
+
+
+def read_griddata(
+    f: IO[str],
+    simroot: Path,
+    sections: Dict[str, Tuple[type, Callable]],
+    shape: Tuple[int],
+) -> Dict[str, dask.array.Array]:
+    """
+    Read GRID DATA section.
+
+    Constants are lazily allocated; external files are lazily read.  Internal
+    arrays are eagerly read, but converted to a dask array for consistency.
+
+    Parameters
+    ----------
+    f: IO[str]
+    simroot: Path
+        Root path of simulation. Used for reading external files.
+    sections: Dict[str, (type, Callabe)]
+        Dictionary containing type of grid data entry, and function to compute
+        number of max_rows to read.
+    shape: Tuple[int]
+        Full shape of model read from DIS file. Should have length:
+            * DIS: 3
+            * DISV: 2
+            * DISU: 1
+
+    Returns
+    -------
+    variables: Dict[str, xr.DataArray]
+    """
+    content = {}
+    for section, (dtype, compute_max_rows) in sections.items():
+        try:
+            layered = advance_to_griddata_section(f, section)
+            max_rows, read_shape = compute_max_rows(shape, layered)
+            kwargs = {
+                "f": f,
+                "simroot": simroot,
+                "dtype": dtype,
+                "shape": read_shape,
+                "max_rows": max_rows,
+            }
+
+            if layered:
+                nlayer = shape[0]
+                data = dask.array.stack(
+                    [read_array(**kwargs) for _ in range(nlayer)],
+                    axis=0,
+                )
+            else:
+                data = read_array(**kwargs)
+
+            content[section] = data
+
+        except Exception as e:
+            raise type(e)(f"{e}\n Error reading {section} in {f.name}") from e
+
+    return content
diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
new file mode 100644
index 000000000..5a59e99bc
--- /dev/null
+++ b/imod/mf6/read_input/list_input.py
@@ -0,0 +1,102 @@
+"""
+Utilities for reading MODFLOW6 list input.
+"""
+from pathlib import Path
+from typing import IO, List, Tuple
+
+import dask.array
+import numpy as np
+
+from .common import (
+    read_external_binaryfile,
+    read_external_textfile,
+    read_internal,
+    strip_line,
+)
+
+
+def recarr_to_dense(recarr, index_columns, fields, shape) -> List[np.ndarray]:
+    """
+    Convert a record array to separate numpy arrays. Uses the index columns to
+    place the values in a dense array form.
+    """
+    # MODFLOW6 is 1-based, Python is 0-based
+    indices = [recarr[column] - 1 for column in index_columns]
+    variables = []
+    for field in fields:
+        data = np.full(shape, np.nan)
+        data[indices] = recarr[field]
+        variables.append(data)
+    return variables
+
+
+def read_internal_listinput(
+    f: IO[str],
+    dtype: type,
+    index_columns: Tuple[str],
+    fields: Tuple[str],
+    shape: Tuple[int],
+    max_rows: int,
+) -> List[dask.array.Array]:
+    recarr = read_internal(f, max_rows, dtype)
+    return recarr_to_dense(recarr, index_columns, fields, shape)
+
+
+def read_external_listinput(
+    path: Path,
+    dtype: type,
+    index_columns: Tuple[str],
+    fields: Tuple[str],
+    shape: Tuple[int],
+    binary: bool,
+):
+    """
+    Read external list input, separate and reshape to a dense array form.
+    """
+    if binary:
+        recarr = read_external_binaryfile(path, dtype)
+    else:
+        recarr = read_external_textfile(path, dtype)
+    return recarr_to_dense(recarr, index_columns, fields, shape)
+
+
+def read_listinput(
+    f: IO[str],
+    simroot: Path,
+    dtype: type,
+    index_columns: Tuple[str],
+    fields: Tuple[str],
+    shape: Tuple[int],
+    max_rows: int,
+) -> List[dask.array.Array]:
+    """
+    MODFLOW6 list input functionality.
+    """
+    # Store position so week can move back in the file if data is stored
+    # internally.
+    position = f.tell()
+
+    # Read and check the first line.
+    firstline = f.readline()
+    stripped = strip_line(firstline)
+    separated = stripped.split()
+    first = separated[0]
+
+    if first == "open/close":
+        fname = separated[1]
+        binary = "(binary)" in stripped
+        path = simroot / fname
+
+        nout = len(fields)
+        x = dask.delayed(read_external_listinput, nout=nout)(
+            path, dtype, index_columns, fields, shape, binary
+        )
+        variable_values = [
+            dask.array.from_delayed(a, shape=shape, dtype=dtype) for a in x
+        ]
+    else:
+        f.seek(position)
+        x = read_internal_listinput(f, dtype, index_columns, fields, shape, max_rows)
+        variable_values = [dask.array.from_array(a) for a in x]
+
+    return variable_values
diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py
index f14b23c08..6e4680b4e 100644
--- a/imod/mf6/simulation.py
+++ b/imod/mf6/simulation.py
@@ -1,6 +1,7 @@
 import collections
 import pathlib
 import subprocess
+from typing import Union
 
 import jinja2
 import numpy as np
@@ -8,6 +9,8 @@
 
 import imod
 
+from .read_input import read_sim_namefile, read_tdis, tdis_time
+
 
 class Modflow6Simulation(collections.UserDict):
     def _initialize_template(self):
@@ -93,6 +96,39 @@ def render(self):
         d["solutiongroups"] = [[("ims6", f"{solvername}.ims", modelnames)]]
         return self._template.render(d)
 
+    @classmethod
+    def open(cls, path: Union[str, pathlib.Path]):
+        simroot = pathlib.Path(path)
+        name = simroot.stem
+        content = read_sim_namefile(path)
+
+        # Get the times from the time discretization file.
+        tdis_path = content["tdis6"]
+        tdis = read_tdis(tdis_path)
+        globaltimes = tdis_time(tdis)
+
+        # Initialize the simulation
+        simulation = cls(name=name, **content)
+
+        # Collect the models
+        MODEL_CLASSES = {
+            "gwf6": imod.mf6.GroundwaterFlowModel,
+        }
+        models = {}
+        for i, (ftype, fname, pname) in enumerate(content["models"]):
+            if pname is None:
+                modelname = f"{ftype}_{i+1}"
+            # Check for duplicate entry
+            elif pname in models:
+                modelname = f"{pname}_{i+1}"
+            else:
+                modelname = pname
+
+            _class = MODEL_CLASSES[ftype]
+            simulation[modelname] = _class.open(fname, simroot, globaltimes)
+
+        return simulation
+
     def write(self, directory=".", binary=True):
         # Check models for required content
         for key, model in self.items():

From e1f2856c39ac4b1136b597d1522fdb6de420d17b Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 4 Apr 2022 20:02:36 +0200
Subject: [PATCH 04/35] Simulation.open() seems to be running.

Introduce StorageBase for Storage packages to avoid repetition
---
 imod/mf6/model.py                 |  77 +++++++--------
 imod/mf6/pkgbase.py               |  59 ++++++++++--
 imod/mf6/read_input/__init__.py   | 152 +++++++++++++++++++++++++++---
 imod/mf6/read_input/grid_data.py  |  90 +++++++++++++-----
 imod/mf6/read_input/list_input.py |  20 +++-
 imod/mf6/simulation.py            |   4 +-
 imod/mf6/sto.py                   | 123 +++++++++++++++---------
 7 files changed, 392 insertions(+), 133 deletions(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 4208fb61d..da457add9 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -31,6 +31,8 @@ class GroundwaterFlowModel(Model):
 
     @staticmethod
     def _PACKAGE_CLASSES() -> Dict[str, Type]:
+        # mf6.OutputControl is skipped
+        # mf6.StorageCoefficient handled by SpecificStorage
         return {
             package._pkg_id: package
             for package in (
@@ -41,16 +43,12 @@ def _PACKAGE_CLASSES() -> Dict[str, Type]:
                 mf6.Evapotranspiration,
                 mf6.GeneralHeadBoundary,
                 mf6.InitialConditions,
-                mf6.Solution,
                 mf6.NodePropertyFlow,
-                mf6.OutputControl,
                 mf6.Recharge,
                 mf6.River,
                 mf6.SpecificStorage,
-                mf6.Storage,
-                mf6.StorageCoefficient,
-                mf6.WellDisStructured,
-                mf6.WellDisVertices,
+                # mf6.WellDisStructured,
+                # mf6.WellDisVertices,
             )
         }
 
@@ -146,7 +144,10 @@ def open(
         globaltimes: np.ndarray,
     ):
         content = read_gwf_namefile(simroot / path)
-        model = cls(**content)
+        model = cls(
+            newton=content.get("newton", False),
+            under_relaxation=content.get("under_relaxation", False),
+        )
 
         # Search for the DIS/DISV/DISU package first. This provides us with
         # the coordinates and dimensions to instantiate the other packages.
@@ -154,56 +155,50 @@ def open(
         dis_packages = [
             tup for tup in content["packages"] if tup[0] in ("dis6", "disv6", "disu6")
         ]
-        packages = [
-            tup
-            for tup in content["packages"]
-            if tup[0] not in ("dis6", "disv6", "disu6")
-        ]
-
-        if len(dis_packages) > 1:
-            raise ValueError(f"More than one DIS/DISV/DISU package in {path}")
+        if len(dis_packages) != 1:
+            raise ValueError(f"Exactly one DIS/DISV/DISU package required in {path}")
 
         disftype, disfname, dispname = dis_packages[0]
         diskey = disftype[:-1]
         package = classes[diskey]
-        if dispname is None:
-            disname = diskey
-
         dis_package = package.open(
             disfname,
             simroot,
         )
-        model[disname] = dis_package
+        model[dispname] = dis_package
         shape = dis_package["idomain"].shape
-        coords = dis_package["idomain"].coords
+        coords = dict(dis_package["idomain"].coords)
         dims = dis_package["idomain"].dims
 
-        # Now read the rest of the packages.
-        seen = set()
-        for i, (ftype, fname, pname) in packages:
-
-            key = ftype[:-1]  # Remove the last number (riv6 -> riv).
-            package = classes[key]
-
-            # Fill in a name if none is given.
-            if pname is None:
-                pname = key
-
-            # Ensure a unique name is generated.
-            if pname in seen:
-                pkgname = f"{pname}_{i+1}"
+        # Non-dis packages
+        packages = [
+            tup
+            for tup in content["packages"]
+            if tup[0] not in ("dis6", "disv6", "disu6")
+        ]
+        # Make sure names are unique by counting first, and appending numbers
+        # if they are not unique.
+        occurence = collections.Counter([tup[2] for tup in packages])
+        counter = collections.defaultdict(int)
+        for ftype, fname, pname in packages:
+            key = ftype[:-1]
+            package = classes.get(key)
+            if package is None:
+                continue
+            if occurence[pname] > 1:
+                pkgname = f"{pname}_{counter[pname]}"
+                counter[pname] += 1
             else:
                 pkgname = pname
-            seen.add(pkgname)
 
             # Create the package and add it to the model.
             model[pkgname] = package.open(
-                fname,
-                simroot,
-                shape,
-                coords,
-                dims,
-                globaltimes,
+                path=fname,
+                simroot=simroot,
+                shape=shape,
+                coords=coords,
+                dims=dims,
+                globaltimes=globaltimes,
             )
 
         return model
diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index ef43b11af..73a4e1a6d 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -1,12 +1,14 @@
 import abc
+import inspect
 import pathlib
+from typing import Any, Dict
 
 import jinja2
 import numpy as np
 import xarray as xr
 import xugrid as xu
 
-from .read_input import get_function_kwargnames, read_blockfile
+from .read_input import read_blockfile, read_boundary_blockfile, shape_to_max_rows
 
 
 def dis_recarr(arrdict, layer, notnull):
@@ -352,6 +354,51 @@ def write(self, directory, pkgname, globaltimes, binary):
                             path = pkgdirectory / f"{key}.dat"
                             self.write_text_griddata(path, da, dtype)
 
+    @classmethod
+    def filter_and_rename(cls, content: Dict[str, Any]) -> Dict[str, Any]:
+        """
+        Filter content for keywords that are required by the ``__init__`` function
+        of a class.
+
+        Rename keywords to the expected form for the class, as given by
+        ``_keyword_map``.
+
+        Arguments
+        ---------
+        init: Callable
+            __init__ of the
+        """
+        inverse_keywords = {v: k for k, v in cls._keyword_map.items()}
+        init_kwargs = inspect.getfullargspec(cls.__init__).args
+
+        filtered = {}
+        for key, value in content.items():
+            newkey = inverse_keywords.get(key, key)
+            if newkey in init_kwargs:
+                filtered[newkey] = value
+
+        return filtered
+
+    @classmethod
+    def open(cls, path, simroot, shape, coords, dims, globaltimes):
+        sections = {
+            cls._keyword_map.get(field, field): (dtype, shape_to_max_rows)
+            for field, dtype in cls._grid_data.items()
+        }
+        content = read_blockfile(
+            simroot / path,
+            simroot=simroot,
+            sections=sections,
+            shape=shape,
+        )
+
+        griddata = content.pop("griddata")
+        for field, data in griddata.items():
+            content[field] = xr.DataArray(data, coords, dims)
+
+        filtered_content = cls.filter_and_rename(content)
+        return cls(**filtered_content)
+
     def _netcdf_path(self, directory, pkgname):
         """create path for netcdf, this function is also used to create paths to use inside the qgis projectfiles"""
         return directory / pkgname / f"{self._pkg_id}.nc"
@@ -529,12 +576,11 @@ def write(self, directory, pkgname, globaltimes, binary):
 
     @classmethod
     def open(cls, path, simroot, shape, coords, dims, globaltimes):
-        content = read_blockfile(
+        content = read_boundary_blockfile(
             simroot / path,
             simroot,
             fields=cls._period_data,
             shape=shape,
-            dims=dims,
         )
 
         period_index = content.pop("period_index")
@@ -546,12 +592,7 @@ def open(cls, path, simroot, shape, coords, dims, globaltimes):
         for field, data in period_data.items():
             content[field] = xr.DataArray(data, coords, dims)
 
-        filtered_content = {
-            k: v
-            for k, v in content.items()
-            if k in get_function_kwargnames(cls.__init__)
-        }
-
+        filtered_content = cls.filter_and_rename(content)
         return cls(**filtered_content)
 
 
diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index cb6e3e3b9..10e6b7ab4 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -63,7 +63,7 @@ def parse_package(stripped: str, fname: str) -> Tuple[str, str, str]:
     nwords = len(separated)
     if nwords == 2:
         ftype, fname = separated
-        pname = ftype
+        pname = ftype[:-1]  # split the number, e.g. riv6 -> riv
     elif nwords == 3:
         ftype, fname, pname = separated
     else:
@@ -129,6 +129,10 @@ def delc_max_rows(shape, layered) -> int:
 
 
 def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
+    """
+    Use start_date, time_units, and period duration to create datetime
+    timestaps for the stress periods.
+    """
     # https://numpy.org/doc/stable/reference/arrays.datetime.html#datetime-units
     TIME_UNITS = {
         "unknown": None,
@@ -138,7 +142,7 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
         "days": "D",
         "years": "Y",
     }
-    unit = TIME_UNITS.get(tdis["units"])
+    unit = TIME_UNITS.get(tdis["time_units"])
 
     start = None
     if "start_date_time" in tdis:
@@ -153,14 +157,14 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
         times = np.full(timedelta.size, start)
         times[1:] += timedelta[:-1]
     else:
-        possibilities = ",".join(list(TIME_UNITS.keys())[1:])
+        possibilities = ", ".join(list(TIME_UNITS.keys())[1:])
         warnings.warn(
-            "Cannot convert time coordinate to datetime. "
-            "Falling back to (unitless) floating point time coordinate. "
+            "Cannot convert time coordinate to datetime."
+            "Falling back to (unitless) floating point time coordinate.\n"
             f"time_units should be one of: {possibilities}; "
             "start_date_time should be according to ISO 8601."
         )
-        times = np.full(timedelta.size, 0.0)
+        times = np.full(cumulative_length.size, 0.0)
         times[1:] += cumulative_length[:-1]
 
     return times
@@ -201,6 +205,44 @@ def read_gwf_namefile(path: Path) -> Dict[str, Any]:
     return content
 
 
+def read_blockfile(
+    path: Path,
+    simroot: Path,
+    sections: Dict[str, Tuple[type, Callable]],
+    shape: Tuple[int],
+) -> Dict[str, Any]:
+    """
+    Read blockfile of a "standard" MODFLOW6 package: NPF, IC, etc.
+
+    External files are lazily read using dask, constants are lazily allocated,
+    and internal values are eagerly read and then converted to dask arrays for
+    consistency.
+
+    Parameters
+    ----------
+    path: Path
+    simroot: Path
+        Root path of the simulation, used for reading external files.
+    sections: Dict[str, Tuple[type, Callable]]
+        Contains for every array entry its type, and a function to compute from
+        shape the number of rows to read in case of internal values.
+    shape: Tuple[int]
+        DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
+
+    Returns
+    -------
+    content: Dict[str, Any]
+        Content of the block file. Grid data arrays are stored as dask arrays.
+    """
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "griddata")
+        content["griddata"] = read_griddata(f, simroot, sections, shape)
+
+    return content
+
+
 def read_package_periods(
     f: IO[str],
     simroot: Path,
@@ -214,11 +256,30 @@ def read_package_periods(
     Read blockfile periods section of a "standard" MODFLOW6 package: RIV, DRN,
     GHB, etc.
 
+    External files are lazily read using dask and internal values are eagerly
+    read and then converted to dask arrays for consistency.
+
     Parameters
     ----------
     f: IO[str]
         File handle.
     simroot:
+        Root path of the simulation, used for reading external files.
+    dtype: dtype
+        Generally a numpy structured dtype.
+    fields: List[str]
+        The fields (generally np.float64) of the dtype.
+    index_columns: List[str]
+        The index columns (np.int32) of the dtype.
+    max_rows: int
+        Number of rows to read in case of internal values.
+    shape: Tuple[int]
+        DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
+
+    Returns
+    -------
+    period_index: np.ndarray of integers
+    variable_values: Dict[str, dask.array.Array]
     """
     periods = []
     dask_lists = defaultdict(list)
@@ -249,20 +310,37 @@ def read_package_periods(
         field: dask.array.stack(dask_list, axis=0)
         for field, dask_list in dask_lists.items()
     }
-    return periods, variable_values
+    return np.array(periods) - 1, variable_values
 
 
-def read_blockfile(
+def read_boundary_blockfile(
     path: Path,
     simroot: Path,
     fields: Tuple[str],
     shape: Tuple[int],
-    dims: Tuple[str],
 ) -> Dict[str, Any]:
     """
-    Read blockfile of a "standard" MODFLOW6 package: RIV, DRN, GHB, etc.
+    Read blockfile of a "standard" MODFLOW6 boundary condition package: RIV,
+    DRN, GHB, etc.
+
+    External files are lazily read using dask and internal values are eagerly
+    read and then converted to dask arrays for consistency.
+
+    Parameters
+    ----------
+    path: Path
+    simroot: Path
+        Root path of the simulation, used for reading external files.
+    fields: List[str]
+        The fields (generally np.float64) of the dtype.
+    shape: Tuple[int]
+        DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
+
+    Returns
+    -------
+    content: Dict[str, Any]
     """
-    ndim = len(dims)
+    ndim = len(shape)
     index_columns = {
         1: ("node",),
         2: ("layer", "cell2d"),
@@ -294,5 +372,53 @@ def read_blockfile(
     return content
 
 
-def get_function_kwargnames(f: Callable) -> List[str]:
-    return inspect.getfullargspec(f).args
+def read_sto_blockfile(
+    path: Path,
+    simroot: Path,
+    sections: Dict[str, Tuple[type, Callable]],
+    shape: Tuple[int],
+) -> Dict[str, Any]:
+    """
+    Read blockfile of MODFLOW6 Storage package.
+
+    Parameters
+    ----------
+    path: Path
+    simroot: Path
+        Root path of the simulation, used for reading external files.
+    sections: Dict[str, Tuple[type, Callable]]
+        Contains for every array entry its type, and a function to compute from
+        shape the number of rows to read in case of internal values.
+    shape: Tuple[int]
+        DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
+
+    Returns
+    -------
+    content: Dict[str, Any]
+        Content of the block file. Grid data arrays are stored as dask arrays.
+    """
+    with open(path, "r") as f:
+        advance_to_header(f, "options")
+        content = read_key_value_block(f, parse_option)
+        advance_to_header(f, "griddata")
+        content["griddata"] = read_griddata(f, simroot, sections, shape)
+
+        periods = {}
+        key = advance_to_period(f)
+        while key is not None:
+            entry = read_key_value_block(f, parse_option)
+            value = next(iter(entry.keys()))
+            if value == "steady-state":
+                value = False
+            elif value == "transient":
+                value = True
+            else:
+                raise ValueError(
+                    f'Expected "steady-state" or "transient". Received: {value}'
+                )
+            periods[key] = value
+            key = advance_to_period(f)
+
+        content["periods"] = periods
+
+    return content
diff --git a/imod/mf6/read_input/grid_data.py b/imod/mf6/read_input/grid_data.py
index 9623d62d6..3afa3a33b 100644
--- a/imod/mf6/read_input/grid_data.py
+++ b/imod/mf6/read_input/grid_data.py
@@ -17,27 +17,54 @@
 )
 
 
-def advance_to_griddata_section(f: IO[str], section: str) -> None:
+def advance_to_griddata_section(f: IO[str]) -> Tuple[str, bool]:
     line = None
-    # Empty line is end-of-file
+    # Empty line is end-of-file.
+    # Should always exit early in while loop, err otherwise.
     while not end_of_file(line):
         line = f.readline()
         stripped = strip_line(line)
-        # Return if start has been found
-        if stripped == section:
-            return False
-        elif stripped == f"{section} layered":
-            return True
-        # Continue if line is empty
-        elif stripped == "":
+        if stripped == "":
             continue
-        # Raise if the line has (unexpected) content
+        elif "end" in stripped:
+            return None, False
+        elif "layered" in stripped:
+            layered = True
+            section = stripped.split()[0]
+            return section, layered
         else:
-            break
-    raise ValueError(f'"{section}" is not present in file {f.name}')
+            layered = False
+            section = stripped
+            return section, layered
+    raise ValueError(f"No end of griddata specified in {f.name}")
 
 
 def shape_to_max_rows(shape: Tuple[int], layered: bool) -> int:
+    """
+    Compute the number of rows to read in case the data is internal. In case
+    of DIS, the number of (table) rows equals the number of layers times the
+    number of (model) rows; a single row then contains ncolumn values.
+
+    A DISV does not have rows and columns, the number of rows is equal to the
+    number of layers.
+
+    In case of LAYERED, each layer is written on a separate row.
+
+    In case of DISU, all values are written one a single row, and LAYERED input
+    is not allowed.
+
+    Parameters
+    ----------
+    shape: Tuple[int]
+    layered: bool
+
+    Parameters
+    ----------
+    max_rows: int
+        Reduced number if layered is True.
+    read_shape: Tuple[int]
+        First dimension removed if layered is True.
+    """
     ndim = len(shape)
     if ndim == 3:
         nlayer, nrow, _ = shape
@@ -99,12 +126,22 @@ def read_array(
     """
     MODFLOW6 READARRAY functionality for grid data.
 
+    External files are lazily read using dask, constants are lazily allocated,
+    and internal values are eagerly read and then converted to dask arrays for
+    consistency.
+
     Parameters
     ----------
+    f: IO[str]
+    simroot: Path
+    dtype: type
+    max_rows: int
+        Number of rows to read in case of internal data values.
+    shape: Tuple[int]
 
     Returns
     -------
-    array: dask.array.Array
+    array: dask.array.Array of size ``shape``
     """
     firstline = f.readline()
     stripped = strip_line(firstline)
@@ -146,8 +183,9 @@ def read_griddata(
     """
     Read GRID DATA section.
 
-    Constants are lazily allocated; external files are lazily read.  Internal
-    arrays are eagerly read, but converted to a dask array for consistency.
+    External files are lazily read using dask, constants are lazily allocated,
+    and internal values are eagerly read and then converted to dask arrays for
+    consistency.
 
     Parameters
     ----------
@@ -165,12 +203,19 @@ def read_griddata(
 
     Returns
     -------
-    variables: Dict[str, xr.DataArray]
+    variables: Dict[str, dask.array.Array]
     """
     content = {}
-    for section, (dtype, compute_max_rows) in sections.items():
-        try:
-            layered = advance_to_griddata_section(f, section)
+    try:
+        section, layered = advance_to_griddata_section(f)
+        while section is not None:
+            try:
+                dtype, compute_max_rows = sections[section]
+            except KeyError:
+                raise ValueError(
+                    f"Unexpected section: {section}. "
+                    f"Expected one of: {', '.join(sections.keys())}"
+                )
             max_rows, read_shape = compute_max_rows(shape, layered)
             kwargs = {
                 "f": f,
@@ -191,7 +236,10 @@ def read_griddata(
 
             content[section] = data
 
-        except Exception as e:
-            raise type(e)(f"{e}\n Error reading {section} in {f.name}") from e
+            # Advance to next section
+            section, layered = advance_to_griddata_section(f)
+
+    except Exception as e:
+        raise type(e)(f"{e}\n Error reading {f.name}") from e
 
     return content
diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
index 5a59e99bc..a7d026eff 100644
--- a/imod/mf6/read_input/list_input.py
+++ b/imod/mf6/read_input/list_input.py
@@ -70,7 +70,25 @@ def read_listinput(
     max_rows: int,
 ) -> List[dask.array.Array]:
     """
-    MODFLOW6 list input functionality.
+    MODFLOW6 list input reading functionality.
+
+    MODFLOW6 list input is "sparse": it consists of a cell id and a number of
+    values. Depending on whether the model is discretized according to DIS,
+    DISV, or DISU; this cell id may be a tuple of size 3, 2, or 1.
+
+    Parameters
+    ----------
+    f: IO[str]
+        File handle.
+    simroot: Path
+        Root path of simulation. Used for reading external files.
+    dtype: type
+
+    index_columns: Tuple[str]
+    fields: Tuple[str]
+    shape: Tuple[int]
+    max_rows: int
+
     """
     # Store position so week can move back in the file if data is stored
     # internally.
diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py
index 6e4680b4e..7cef7a42b 100644
--- a/imod/mf6/simulation.py
+++ b/imod/mf6/simulation.py
@@ -98,7 +98,7 @@ def render(self):
 
     @classmethod
     def open(cls, path: Union[str, pathlib.Path]):
-        simroot = pathlib.Path(path)
+        simroot = pathlib.Path(path).parent
         name = simroot.stem
         content = read_sim_namefile(path)
 
@@ -108,7 +108,7 @@ def open(cls, path: Union[str, pathlib.Path]):
         globaltimes = tdis_time(tdis)
 
         # Initialize the simulation
-        simulation = cls(name=name, **content)
+        simulation = cls(name=name)
 
         # Collect the models
         MODEL_CLASSES = {
diff --git a/imod/mf6/sto.py b/imod/mf6/sto.py
index c5e7593d8..a5728ee05 100644
--- a/imod/mf6/sto.py
+++ b/imod/mf6/sto.py
@@ -1,7 +1,13 @@
+from pathlib import Path
+from typing import Tuple
+
 import numpy as np
+import xarray as xr
 
 from imod.mf6.pkgbase import Package
 
+from .read_input import read_sto_blockfile, shape_to_max_rows
+
 
 class Storage(Package):
     def __init__(*args, **kwargs):
@@ -10,7 +16,72 @@ def __init__(*args, **kwargs):
         )
 
 
-class SpecificStorage(Package):
+class StorageBase(Package):
+    def _render(self, directory, pkgname, globaltimes, binary):
+        d = {}
+        stodirectory = directory / "sto"
+        for varname in self._grid_data.keys():
+            key = self._keyword_map.get(varname, varname)
+            layered, value = self._compose_values(
+                self[varname], stodirectory, key, binary=binary
+            )
+            if self._valid(value):  # skip False or None
+                d[f"{key}_layered"], d[key] = layered, value
+
+        periods = {}
+        if "time" in self.dataset["transient"].coords:
+            package_times = self.dataset["transient"].coords["time"].values
+            starts = np.searchsorted(globaltimes, package_times) + 1
+            for i, s in enumerate(starts):
+                periods[s] = self.dataset["transient"].isel(time=i).values[()]
+        else:
+            periods[1] = self.dataset["transient"].values[()]
+
+        d["periods"] = periods
+        return d
+
+    @staticmethod
+    def open(
+        path: Path,
+        simroot: Path,
+        shape: Tuple[int],
+        coords: Tuple[str],
+        dims: Tuple[str],
+        globaltimes: np.ndarray,
+    ):
+        sections = {
+            "iconvert": (np.int32, shape_to_max_rows),
+            "ss": (np.float64, shape_to_max_rows),
+            "sy": (np.float64, shape_to_max_rows),
+        }
+        content = read_sto_blockfile(
+            path=path,
+            simroot=simroot,
+            sections=sections,
+            shape=shape,
+        )
+
+        griddata = content.pop("griddata")
+        for field, data in griddata.items():
+            content[field] = xr.DataArray(data, coords, dims)
+        periods = content.pop("periods")
+        content["transient"] = xr.DataArray(
+            list(periods.values()),
+            coords={"time": globaltimes[list(periods.keys())]},
+            dims=("time",),
+        )
+
+        storagecoefficient = content.pop("storagecoefficient", False)
+        if storagecoefficient:
+            cls = StorageCoefficient
+        else:
+            cls = SpecificStorage
+
+        filtered_content = cls.filter_and_rename(content)
+        return cls(**filtered_content)
+
+
+class SpecificStorage(StorageBase):
     """
     Storage Package with specific storage.
 
@@ -64,34 +135,14 @@ def __init__(self, specific_storage, specific_yield, transient, convertible):
         self.dataset["transient"] = transient
 
     def render(self, directory, pkgname, globaltimes, binary):
-        d = {}
-        stodirectory = directory / "sto"
-        for varname in ["specific_storage", "specific_yield", "convertible"]:
-            key = self._keyword_map.get(varname, varname)
-            layered, value = self._compose_values(
-                self[varname], stodirectory, key, binary=binary
-            )
-            if self._valid(value):  # skip False or None
-                d[f"{key}_layered"], d[key] = layered, value
-
-        periods = {}
-        if "time" in self.dataset["transient"].coords:
-            package_times = self.dataset["transient"].coords["time"].values
-            starts = np.searchsorted(globaltimes, package_times) + 1
-            for i, s in enumerate(starts):
-                periods[s] = self.dataset["transient"].isel(time=i).values[()]
-        else:
-            periods[1] = self.dataset["transient"].values[()]
-
-        d["periods"] = periods
-
+        d = self._render(directory, globaltimes, binary)
         return self._template.render(d)
 
 
-class StorageCoefficient(Package):
+class StorageCoefficient(StorageBase):
     """
-    Storage Package with a storage coefficient.  Be careful,
-    this is not the same as the specific storage.
+    Storage Package with a storage coefficient. Be careful, this is not the
+    same as the specific storage.
 
     From wikipedia (https://en.wikipedia.org/wiki/Specific_storage):
 
@@ -153,26 +204,6 @@ def __init__(self, storage_coefficient, specific_yield, transient, convertible):
         self.dataset["transient"] = transient
 
     def render(self, directory, pkgname, globaltimes, binary):
-        d = {}
-        stodirectory = directory / "sto"
-        for varname in ["storage_coefficient", "specific_yield", "convertible"]:
-            key = self._keyword_map.get(varname, varname)
-            layered, value = self._compose_values(
-                self[varname], stodirectory, key, binary=binary
-            )
-            if self._valid(value):  # skip False or None
-                d[f"{key}_layered"], d[key] = layered, value
-
-        periods = {}
-        if "time" in self.dataset["transient"].coords:
-            package_times = self.dataset["transient"].coords["time"].values
-            starts = np.searchsorted(globaltimes, package_times) + 1
-            for i, s in enumerate(starts):
-                periods[s] = self.dataset["transient"].isel(time=i).values[()]
-        else:
-            periods[1] = self.dataset["transient"].values[()]
-
-        d["periods"] = periods
+        d = self._render(directory, globaltimes, binary)
         d["storagecoefficient"] = True
-
         return self._template.render(d)

From c1d1e7583d39dc10034541e31880c288b7137a4d Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 6 Apr 2022 21:22:19 +0200
Subject: [PATCH 05/35] Adjustments:

* Use pandas for ~3 times faster reading of list input
* Add logic for mimicking parts of Fortran read intrinsic: double formats (to_float), comma and whitespace separation (split_string), weird concatenation operators in text files (read_fortran_deflated_text_array).
* Allow parse_option to return more than one option (required for auxiliary).
* Add auxiliary and boundnames entries to list input dtype.
* Fix delr, delc (both 1D) and top (2D) in read_dis_blockfile.
* Intercept errors and attach file name.
---
 imod/mf6/dis.py                   |  19 +++---
 imod/mf6/model.py                 |  21 ++++---
 imod/mf6/read_input/__init__.py   |  99 ++++++++++++++++++++---------
 imod/mf6/read_input/common.py     | 101 +++++++++++++++++++++++++-----
 imod/mf6/read_input/grid_data.py  |  21 ++++---
 imod/mf6/read_input/list_input.py |  68 ++++++++++++++------
 6 files changed, 240 insertions(+), 89 deletions(-)

diff --git a/imod/mf6/dis.py b/imod/mf6/dis.py
index 0998a9427..8c7a1a98b 100644
--- a/imod/mf6/dis.py
+++ b/imod/mf6/dis.py
@@ -68,15 +68,20 @@ def open(cls, path: str, simroot: Path) -> "StructuredDiscretization":
         ymin = float(content.get("yorigin", 0.0))
         dx = griddata.pop("delr")
         dy = griddata.pop("delc")
-        coords = {
-            "layer": np.arange(1, nlay + 1),
-            "y": ((ymin + np.cumsum(dy)) - 0.5 * dx)[::-1],
-            "x": (xmin + np.cumsum(dx)) - 0.5 * dx,
-        }
-        dims = ("layer", "y", "x")
+        x = (xmin + np.cumsum(dx)) - 0.5 * dx
+        y = ((ymin + np.cumsum(dy)) - 0.5 * dy)[::-1]
+
+        # Top is a 2D array unlike the others
+        top = xr.DataArray(
+            data=griddata.pop("top"),
+            coords={"y": y, "x": x},
+            dims=("y", "x"),
+        )
 
+        coords = {"layer": np.arange(1, nlay + 1), "y": y, "x": x}
+        dims = ("layer", "y", "x")
         inverse_keywords = {v: k for k, v in cls._keyword_map.items()}
-        variables = {}
+        variables = {"top": top}
         for key, value in griddata.items():
             invkey = inverse_keywords.get(key, key)
             variables[invkey] = xr.DataArray(
diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index da457add9..096ec2671 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -191,15 +191,18 @@ def open(
             else:
                 pkgname = pname
 
-            # Create the package and add it to the model.
-            model[pkgname] = package.open(
-                path=fname,
-                simroot=simroot,
-                shape=shape,
-                coords=coords,
-                dims=dims,
-                globaltimes=globaltimes,
-            )
+            try:
+                # Create the package and add it to the model.
+                model[pkgname] = package.open(
+                    path=fname,
+                    simroot=simroot,
+                    shape=shape,
+                    coords=coords,
+                    dims=dims,
+                    globaltimes=globaltimes,
+                )
+            except Exception as e:
+                raise type(e)(f"{e}\n Error reading {fname}") from e
 
         return model
 
diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index 10e6b7ab4..8e9df392c 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -1,11 +1,10 @@
 """
 Utilities for reading MODFLOW6 input files.
 """
-import inspect
 import warnings
 from collections import defaultdict
 from pathlib import Path
-from typing import IO, Any, Callable, Dict, List, Tuple
+from typing import IO, Any, Callable, Dict, List, Tuple, Union
 
 import dask.array
 import numpy as np
@@ -17,13 +16,16 @@
     parse_option,
     read_iterable_block,
     read_key_value_block,
+    split_line,
+    to_float,
 )
 from .grid_data import read_griddata, shape_to_max_rows
 from .list_input import read_listinput
 
 
 def parse_model(stripped: str, fname: str) -> Tuple[str, str, str]:
-    separated = stripped.split()
+    """Parse model entry in the simulation name file."""
+    separated = split_line(stripped)
     nwords = len(separated)
     if nwords == 3:
         return separated
@@ -35,7 +37,8 @@ def parse_model(stripped: str, fname: str) -> Tuple[str, str, str]:
 
 
 def parse_exchange(stripped: str, fname: str) -> Tuple[str, str, str, str]:
-    separated = stripped.split()
+    """Parse exchange entry in the simulation name file."""
+    separated = split_line(stripped)
     nwords = len(separated)
     if nwords == 4:
         return separated
@@ -47,7 +50,11 @@ def parse_exchange(stripped: str, fname: str) -> Tuple[str, str, str, str]:
 
 
 def parse_solutiongroup(stripped: str, fname: str) -> Tuple[str, str]:
-    separated = stripped.split()
+    """Parse solution group entry iyn the simulation name file."""
+    separated = split_line(stripped)
+    if "mxiter" in stripped:
+        return separated
+
     nwords = len(separated)
     if nwords > 2:
         return separated[0], separated[1], separated[2:]
@@ -59,7 +66,8 @@ def parse_solutiongroup(stripped: str, fname: str) -> Tuple[str, str]:
 
 
 def parse_package(stripped: str, fname: str) -> Tuple[str, str, str]:
-    separated = stripped.split()
+    """Parse package entry in model name file."""
+    separated = split_line(stripped)
     nwords = len(separated)
     if nwords == 2:
         ftype, fname = separated
@@ -75,10 +83,11 @@ def parse_package(stripped: str, fname: str) -> Tuple[str, str, str]:
 
 
 def parse_tdis_perioddata(stripped: str, fname: str) -> Tuple[float, int, float]:
-    separated = stripped.split()
+    """Parse a single period data entry in the time discretization file."""
+    separated = split_line(stripped)
     nwords = len(separated)
-    if nwords == 3:
-        return float(separated[0]), int(separated[1]), float(separated[2])
+    if nwords >= 3:
+        return to_float(separated[0]), int(separated[1]), to_float(separated[2])
     else:
         raise ValueError(
             "perlen, nstp, tsmult expected. Received instead: "
@@ -87,6 +96,7 @@ def parse_tdis_perioddata(stripped: str, fname: str) -> Tuple[float, int, float]
 
 
 def read_tdis(path: Path) -> Dict[str, Any]:
+    """Read and parse the content of the time discretization file."""
     with open(path, "r") as f:
         advance_to_header(f, "options")
         content = read_key_value_block(f, parse_option)
@@ -98,20 +108,30 @@ def read_tdis(path: Path) -> Dict[str, Any]:
 
 
 def read_dis_blockfile(path: Path, simroot: Path) -> Dict[str, Any]:
-    def delr_max_rows(shape, layered) -> Tuple[int, int]:
+    """Read and parse structured discretization file."""
+    ROW = 1
+    COLUMN = 2
+
+    def delr_max_rows(shape, layered) -> Tuple[int, Tuple[int]]:
         if layered:
             raise ValueError(f"DELR section in {path} is LAYERED")
-        return 1, shape[2]
+        return 1, (shape[COLUMN],)
 
-    def delc_max_rows(shape, layered) -> int:
+    def delc_max_rows(shape, layered) -> Tuple[int, Tuple[int]]:
         if layered:
             raise ValueError(f"DELC section in {path} is LAYERED")
-        return 1, shape[1]
+        return 1, (shape[ROW],)
+
+    def top_max_rows(shape, layered) -> Tuple[int, Tuple[int]]:
+        if layered:
+            raise ValueError(f"TOP section in {path} is LAYERED")
+        _, nrow, ncolumn = shape
+        return nrow, (nrow, ncolumn)
 
     sections = {
         "delr": (np.float64, delr_max_rows),
         "delc": (np.float64, delc_max_rows),
-        "top": (np.float64, shape_to_max_rows),
+        "top": (np.float64, top_max_rows),
         "botm": (np.float64, shape_to_max_rows),
         "idomain": (np.int32, shape_to_max_rows),
     }
@@ -171,6 +191,7 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
 
 
 def read_solver(path: Path) -> Dict[str, str]:
+    """Read and parse content of solver config file (IMS)."""
     with open(path, "r") as f:
         advance_to_header(f, "options")
         options = read_key_value_block(f, parse_option)
@@ -182,6 +203,7 @@ def read_solver(path: Path) -> Dict[str, str]:
 
 
 def read_sim_namefile(path: Path) -> Dict[str, str]:
+    """Read and parse content of simulation name file."""
     with open(path, "r") as f:
         advance_to_header(f, "options")
         content = read_key_value_block(f, parse_option)
@@ -197,6 +219,7 @@ def read_sim_namefile(path: Path) -> Dict[str, str]:
 
 
 def read_gwf_namefile(path: Path) -> Dict[str, Any]:
+    """Read and parse content of groundwater flow name file."""
     with open(path, "r") as f:
         advance_to_header(f, "options")
         content = read_key_value_block(f, parse_option)
@@ -246,15 +269,15 @@ def read_blockfile(
 def read_package_periods(
     f: IO[str],
     simroot: Path,
-    dtype: type,
-    fields: Tuple[str],
+    dtype: Union[type, np.dtype],
     index_columns: Tuple[str],
+    fields: Tuple[str],
     max_rows: int,
     shape: Tuple[int],
 ) -> Tuple[List[int], Dict[str, dask.array.Array]]:
     """
-    Read blockfile periods section of a "standard" MODFLOW6 package: RIV, DRN,
-    GHB, etc.
+    Read blockfile periods section of a "standard" MODFLOW6 boundary
+    conditions: RIV, DRN, GHB, etc.
 
     External files are lazily read using dask and internal values are eagerly
     read and then converted to dask arrays for consistency.
@@ -265,12 +288,12 @@ def read_package_periods(
         File handle.
     simroot:
         Root path of the simulation, used for reading external files.
-    dtype: dtype
+    dtype: Union[type, np.dtype]
         Generally a numpy structured dtype.
-    fields: List[str]
-        The fields (generally np.float64) of the dtype.
-    index_columns: List[str]
+    index_columns: Tuple[str]
         The index columns (np.int32) of the dtype.
+    fields: Tuple[str]
+        The fields (generally np.float64) of the dtype.
     max_rows: int
         Number of rows to read in case of internal values.
     shape: Tuple[int]
@@ -331,7 +354,7 @@ def read_boundary_blockfile(
     path: Path
     simroot: Path
         Root path of the simulation, used for reading external files.
-    fields: List[str]
+    fields: Tuple[str]
         The fields (generally np.float64) of the dtype.
     shape: Tuple[int]
         DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
@@ -348,23 +371,41 @@ def read_boundary_blockfile(
     }.get(ndim)
     if index_columns is None:
         raise ValueError(f"Length of dimensions should be 1, 2, or 3. Received: {ndim}")
-
-    dtype = np.dtype(
-        [(column, np.int32) for column in index_columns]
-        + [(field, np.float64) for field in fields]
-    )
+    dtype_fields = [(column, np.int32) for column in index_columns] + [
+        (field, np.float64) for field in fields
+    ]
 
     with open(path, "r") as f:
         advance_to_header(f, "options")
         content = read_key_value_block(f, parse_option)
+
+        # Create the dtype from the required variables, and with potential
+        # auxiliary variables.
+        auxiliary = content.pop("auxiliary", None)
+        if auxiliary:
+            # Make sure it's iterable in case of a single value
+            if isinstance(auxiliary, str):
+                auxiliary = (auxiliary,)
+            for aux in auxiliary:
+                dtype_fields.append((aux, np.float64))
+
+        boundnames = content.pop("boundnames", False)
+        if boundnames:
+            dtype_fields.append(("boundnames", str))
+
+        # Create np.dtype, make fields and columns immutable.
+        index_columns = tuple(index_columns)
+        fields = tuple(fields)
+        dtype = np.dtype(dtype_fields)
+
         advance_to_header(f, "dimensions")
         dimensions = read_key_value_block(f, parse_dimension)
         content["period_index"], content["period_data"] = read_package_periods(
             f=f,
             simroot=simroot,
             dtype=dtype,
-            fields=fields,
             index_columns=index_columns,
+            fields=fields,
             max_rows=dimensions["maxbound"],
             shape=shape,
         )
diff --git a/imod/mf6/read_input/common.py b/imod/mf6/read_input/common.py
index 24d8b1cfb..c6d11a01d 100644
--- a/imod/mf6/read_input/common.py
+++ b/imod/mf6/read_input/common.py
@@ -17,15 +17,37 @@ def strip_line(line: str) -> str:
     return before.strip().lower()
 
 
-def find_entry(line: str, entry: str, type: type) -> str:
+def flatten(iterable: Any) -> List[Any]:
+    out = []
+    for part in iterable:
+        out.extend(part)
+    return out
+
+
+def split_line(line: str) -> List[str]:
+    # Maybe check: https://stackoverflow.com/questions/36165050/python-equivalent-of-fortran-list-directed-input
+    # Split on comma and whitespace, like a FORTRAN read would do.
+    return flatten([part.split(",") for part in line.split()])
+
+
+def to_float(string: str) -> float:
+    # Fortran float may contain d (e.g. 1.0d0), Python only accepts e.
+    string = string.replace("d", "e")
+    # Fortran may specify exponents without a letter, e.g. 1.0+1 for 1.0e+1
+    if "e" not in string:
+        string = string[0] + string.replace("+", "e+").replace("-", "e-")
+    return float(string)
+
+
+def find_entry(line: str, entry: str, cast: Callable) -> str:
     if entry not in line:
         return None
     else:
         _, after = line.split(entry)
-        return type(after.split()[0])
+        return cast(after.split()[0])
 
 
-def read_internal(f: IO[str], max_rows: int, dtype: type) -> np.ndarray:
+def read_internal(f: IO[str], dtype: type, max_rows: int) -> np.ndarray:
     return np.loadtxt(
         fname=f,
         dtype=dtype,
@@ -33,19 +55,66 @@ def read_internal(f: IO[str], max_rows: int, dtype: type) -> np.ndarray:
     )
 
 
-def read_external_binaryfile(path: Path, dtype: type) -> np.ndarray:
+def read_external_binaryfile(path: Path, dtype: type, max_rows: int) -> np.ndarray:
     return np.fromfile(
         file=path,
         dtype=dtype,
+        count=max_rows,
         sep="",
     )
 
 
-def read_external_textfile(path: Path, dtype: type) -> np.ndarray:
-    return np.loadtxt(
-        fname=path,
-        dtype=dtype,
-    )
+def read_fortran_deflated_text_array(
+    path: Path, dtype: type, max_rows: int
+) -> np.ndarray:
+    """
+    The Fortran read intrinsic is capable of parsing arrays in many forms.
+    One of those is:
+
+        1.0
+        2*2.0
+        3.0
+
+    Which should be interpreted as: [1.0, 2.0, 2.0, 3.0]
+
+    This function attempts this.
+    """
+    out = np.empty(max_rows, dtype)
+    with open(path, "r") as f:
+        lines = [line.strip() for line in f]
+
+    iterable_lines = iter(lines)
+    start = 0
+    while start < max_rows:
+
+        line = next(iterable_lines)
+        if "*" in line:
+            n, value = line.split("*")
+            n = int(n)
+            value = dtype(value)
+        else:
+            n = 1
+            value = dtype(line)
+
+        end = start + n
+        out[start:end] = value
+        start = end
+
+    return out
+
+
+def read_external_textfile(path: Path, dtype: type, max_rows: int) -> np.ndarray:
+    try:
+        return np.loadtxt(
+            fname=path,
+            dtype=dtype,
+            max_rows=max_rows,
+        )
+    except ValueError as e:
+        if str(e).startswith("could not convert string to float"):
+            return read_fortran_deflated_text_array(path, dtype, max_rows)
+        else:
+            raise e
 
 
 def advance_to_header(f: IO[str], section) -> None:
@@ -54,7 +123,7 @@ def advance_to_header(f: IO[str], section) -> None:
     # Empty line is end-of-file
     while not end_of_file(line):
         line = f.readline()
-        stripped = line.strip().lower()
+        stripped = strip_line(line)
         # Return if start has been found
         if stripped == start:
             return
@@ -111,19 +180,19 @@ def read_iterable_block(f: IO[str], parse: Callable) -> List[Any]:
     raise ValueError(f'"end" of block is not present in file {fname}')
 
 
-def parse_option(stripped: str, fname: str) -> Tuple[str, str]:
+def parse_option(stripped: str, fname: str) -> Tuple[str, Any]:
     separated = stripped.split()
     nwords = len(separated)
-    if nwords == 1:
+    if nwords == 0:
+        raise ValueError(f"Cannot parse {stripped} in {fname}")
+    elif nwords == 1:
         key = separated[0]
         value = True
     elif nwords == 2:
         key, value = separated
     else:
-        raise ValueError(
-            "More than two words found in block:"
-            f"{','.join(separated)} in file {fname}"
-        )
+        key = separated[0]
+        value = separated[1:]
     return key, value
 
 
diff --git a/imod/mf6/read_input/grid_data.py b/imod/mf6/read_input/grid_data.py
index 3afa3a33b..4a0780b70 100644
--- a/imod/mf6/read_input/grid_data.py
+++ b/imod/mf6/read_input/grid_data.py
@@ -13,7 +13,9 @@
     read_external_binaryfile,
     read_external_textfile,
     read_internal,
+    split_line,
     strip_line,
+    to_float,
 )
 
 
@@ -30,7 +32,7 @@ def advance_to_griddata_section(f: IO[str]) -> Tuple[str, bool]:
             return None, False
         elif "layered" in stripped:
             layered = True
-            section = stripped.split()[0]
+            section = split_line(stripped)[0]
             return section, layered
         else:
             layered = False
@@ -39,7 +41,7 @@ def advance_to_griddata_section(f: IO[str]) -> Tuple[str, bool]:
     raise ValueError(f"No end of griddata specified in {f.name}")
 
 
-def shape_to_max_rows(shape: Tuple[int], layered: bool) -> int:
+def shape_to_max_rows(shape: Tuple[int], layered: bool) -> Tuple[int, Tuple[int]]:
     """
     Compute the number of rows to read in case the data is internal. In case
     of DIS, the number of (table) rows equals the number of layers times the
@@ -58,7 +60,7 @@ def shape_to_max_rows(shape: Tuple[int], layered: bool) -> int:
     shape: Tuple[int]
     layered: bool
 
-    Parameters
+    Returns
     ----------
     max_rows: int
         Reduced number if layered is True.
@@ -103,16 +105,17 @@ def constant(value: Any, shape: Tuple[int], dtype: type) -> dask.array.Array:
 def read_internal_griddata(
     f: IO[str], dtype: type, shape: Tuple[int], max_rows: int
 ) -> np.ndarray:
-    return read_internal(f, max_rows, dtype).reshape(shape)
+    return read_internal(f=f, dtype=dtype, max_rows=max_rows).reshape(shape)
 
 
 def read_external_griddata(
     path: Path, dtype: type, shape: Tuple[int], binary: bool
 ) -> np.ndarray:
+    max_rows = np.product(shape)
     if binary:
-        a = read_external_binaryfile(path, dtype)
+        a = read_external_binaryfile(path, dtype, max_rows)
     else:
-        a = read_external_textfile(path, dtype)
+        a = read_external_textfile(path, dtype, max_rows)
     return a.reshape(shape)
 
 
@@ -145,18 +148,18 @@ def read_array(
     """
     firstline = f.readline()
     stripped = strip_line(firstline)
-    separated = stripped.split()
+    separated = split_line(stripped)
     first = separated[0]
 
     if first == "constant":
         factor = None
         array = constant(separated[1], shape, dtype)
     elif first == "internal":
-        factor = find_entry(stripped, "factor", float)
+        factor = find_entry(stripped, "factor", to_float)
         a = read_internal_griddata(f, dtype, shape, max_rows)
         array = dask.array.from_array(a)
     elif first == "open/close":
-        factor = find_entry(stripped, "factor", float)
+        factor = find_entry(stripped, "factor", to_float)
         fname = separated[1]
         binary = "(binary)" in stripped
         path = simroot / fname
diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
index a7d026eff..e75a97b07 100644
--- a/imod/mf6/read_input/list_input.py
+++ b/imod/mf6/read_input/list_input.py
@@ -6,16 +6,17 @@
 
 import dask.array
 import numpy as np
+import pandas as pd
 
-from .common import (
-    read_external_binaryfile,
-    read_external_textfile,
-    read_internal,
-    strip_line,
-)
+from .common import read_external_binaryfile, split_line, strip_line
 
 
-def recarr_to_dense(recarr, index_columns, fields, shape) -> List[np.ndarray]:
+def recarr_to_dense(
+    recarr: np.ndarray,
+    index_columns: Tuple[str],
+    fields: Tuple[str],
+    shape: Tuple[int],
+) -> List[np.ndarray]:
     """
     Convert a record array to separate numpy arrays. Uses the index columns to
     place the values in a dense array form.
@@ -30,33 +31,55 @@ def recarr_to_dense(recarr, index_columns, fields, shape) -> List[np.ndarray]:
     return variables
 
 
+def read_text_listinput(
+    path: Path,
+    dtype: np.dtype,
+    max_rows: int,
+) -> np.ndarray:
+    # This function is over three times faster than:
+    # recarr = np.loadtxt(path, dtype, max_rows=max_rows)
+    d = {key: value[0] for key, value in dtype.fields.items()}
+    df = pd.read_csv(
+        path,
+        header=None,
+        dtype=d,
+        names=d.keys(),
+        delim_whitespace=True,
+        comment="#",
+        nrows=max_rows,
+    )
+    return df.to_records(index=False)
+
+
 def read_internal_listinput(
     f: IO[str],
-    dtype: type,
+    dtype: np.dtype,
     index_columns: Tuple[str],
     fields: Tuple[str],
     shape: Tuple[int],
     max_rows: int,
-) -> List[dask.array.Array]:
-    recarr = read_internal(f, max_rows, dtype)
+) -> List[np.ndarray]:
+    # recarr = read_internal(f, max_rows, dtype)
+    recarr = read_text_listinput(f, dtype, max_rows)
     return recarr_to_dense(recarr, index_columns, fields, shape)
 
 
 def read_external_listinput(
     path: Path,
-    dtype: type,
+    dtype: np.dtype,
     index_columns: Tuple[str],
     fields: Tuple[str],
     shape: Tuple[int],
     binary: bool,
-):
+    max_rows: int,
+) -> List[np.ndarray]:
     """
     Read external list input, separate and reshape to a dense array form.
     """
     if binary:
-        recarr = read_external_binaryfile(path, dtype)
+        recarr = read_external_binaryfile(path, dtype, max_rows)
     else:
-        recarr = read_external_textfile(path, dtype)
+        recarr = read_text_listinput(path, dtype, max_rows)
     return recarr_to_dense(recarr, index_columns, fields, shape)
 
 
@@ -82,13 +105,16 @@ def read_listinput(
         File handle.
     simroot: Path
         Root path of simulation. Used for reading external files.
-    dtype: type
-
+    dtype: np.dtype
     index_columns: Tuple[str]
     fields: Tuple[str]
     shape: Tuple[int]
     max_rows: int
 
+    Returns
+    -------
+    variable_values: List[dask.array.Array]
+        A dask array for every entry in ``fields``.
     """
     # Store position so week can move back in the file if data is stored
     # internally.
@@ -97,22 +123,26 @@ def read_listinput(
     # Read and check the first line.
     firstline = f.readline()
     stripped = strip_line(firstline)
-    separated = stripped.split()
+    separated = split_line(stripped)
     first = separated[0]
 
     if first == "open/close":
         fname = separated[1]
-        binary = "(binary)" in stripped
         path = simroot / fname
+        binary = "(binary)" in stripped
+
+        if binary and "boundnames" in dtype.fields:
+            raise ValueError("(BINARY) does not support BOUNDNAMES")
 
         nout = len(fields)
         x = dask.delayed(read_external_listinput, nout=nout)(
-            path, dtype, index_columns, fields, shape, binary
+            path, dtype, index_columns, fields, shape, binary, max_rows
         )
         variable_values = [
             dask.array.from_delayed(a, shape=shape, dtype=dtype) for a in x
         ]
     else:
+        # Move file position back one line, and try reading internal values.
         f.seek(position)
         x = read_internal_listinput(f, dtype, index_columns, fields, shape, max_rows)
         variable_values = [dask.array.from_array(a) for a in x]

From e24d517fb84892683152a98b38cf6849b8cb0c20 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 8 Apr 2022 12:12:20 +0200
Subject: [PATCH 06/35] Use package name for storage directory

---
 imod/mf6/sto.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/sto.py b/imod/mf6/sto.py
index a5728ee05..6f62866e1 100644
--- a/imod/mf6/sto.py
+++ b/imod/mf6/sto.py
@@ -19,7 +19,7 @@ def __init__(*args, **kwargs):
 class StorageBase(Package):
     def _render(self, directory, pkgname, globaltimes, binary):
         d = {}
-        stodirectory = directory / "sto"
+        stodirectory = directory / pkgname
         for varname in self._grid_data.keys():
             key = self._keyword_map.get(varname, varname)
             layered, value = self._compose_values(
@@ -135,7 +135,7 @@ def __init__(self, specific_storage, specific_yield, transient, convertible):
         self.dataset["transient"] = transient
 
     def render(self, directory, pkgname, globaltimes, binary):
-        d = self._render(directory, globaltimes, binary)
+        d = self._render(directory, pkgname, globaltimes, binary)
         return self._template.render(d)
 
 
@@ -204,6 +204,6 @@ def __init__(self, storage_coefficient, specific_yield, transient, convertible):
         self.dataset["transient"] = transient
 
     def render(self, directory, pkgname, globaltimes, binary):
-        d = self._render(directory, globaltimes, binary)
+        d = self._render(directory, pkgname, globaltimes, binary)
         d["storagecoefficient"] = True
         return self._template.render(d)

From 18697e2d84bb575d1dc5a758c870c046c2e991de Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 8 Apr 2022 12:12:38 +0200
Subject: [PATCH 07/35] Intercept NaNs as xarray yreplaces None by NaN when
 writing to file

---
 imod/mf6/pkgbase.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index 3f6b8716d..0960066d4 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -166,13 +166,18 @@ def sel(self):
 
     def _valid(self, value):
         """
-        Filters values that are None, False, or a numpy.bool_ False.
+        Filters values that are None, False, np.nan, or a numpy.bool_ False.
         Needs to be this specific, since 0.0 and 0 are valid values, but are
         equal to a boolean False.
+
+        Intercepting single NaNs is practical because xarray replaces None by
+        NaN when writing.
         """
         # Test singletons
         if value is False or value is None:
             return False
+        elif isinstance(value, float) and np.isnan(value):
+            return False
         # Test numpy bool (not singleton)
         elif isinstance(value, np.bool_) and not value:
             return False
@@ -226,6 +231,7 @@ def to_sparse(self, arrdict, layer):
                 recarr = disu_recarr(arrdict, layer, notnull)
             else:
                 recarr = dis_recarr(arrdict, layer, notnull)
+
         elif isinstance(self.dataset, xu.UgridDataset):
             recarr = disv_recarr(arrdict, layer, notnull)
         else:
@@ -465,7 +471,7 @@ def write_netcdf(self, directory, pkgname, aggregate_layers=False):
 
         spatial_ds.to_netcdf(path)
         return has_dims
-        
+
     def _to_disu(self, numbers):
         structured = self.dataset.expand_dims("layer")
         # Stack will automatically broadcast to layer
@@ -476,7 +482,7 @@ def _to_disu(self, numbers):
         index = np.add.outer(offset, np.arange(ncell_per_layer))
         dataset = dataset.assign_coords(node=numbers[index])
         return self.__class__(**dataset)
- 
+
 
 class BoundaryCondition(Package, abc.ABC):
     """

From c5d1ff4c723b06d25d92e1fad145be0fba98f405 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 6 May 2022 11:42:24 +0200
Subject: [PATCH 08/35] Do not duplicate entry when replacing exponents when
 casting to float

---
 imod/mf6/read_input/common.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/read_input/common.py b/imod/mf6/read_input/common.py
index c6d11a01d..1142b5149 100644
--- a/imod/mf6/read_input/common.py
+++ b/imod/mf6/read_input/common.py
@@ -35,7 +35,7 @@ def to_float(string: str) -> float:
     string = string.replace("d", "e")
     # Fortran may specify exponents without a letter, e.g. 1.0+1 for 1.0e+1
     if "e" not in string:
-        string = string[0] + string.replace("+", "e+").replace("-", "e-")
+        string = string.replace("+", "e+").replace("-", "e-")
     return float(string)
 
 

From 7ae6e09a61598d57a71a582202b7873a8c7a4ada Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:54:09 +0200
Subject: [PATCH 09/35] Write low level DISU fixes

---
 imod/mf6/__init__.py           |   3 +-
 imod/mf6/disu.py               | 167 +++++++++++++++++++--------------
 imod/mf6/pkgbase.py            |  61 +++++++++---
 imod/templates/mf6/gwf-disu.j2 |   3 +
 4 files changed, 151 insertions(+), 83 deletions(-)

diff --git a/imod/mf6/__init__.py b/imod/mf6/__init__.py
index b183fc772..3785b10ca 100644
--- a/imod/mf6/__init__.py
+++ b/imod/mf6/__init__.py
@@ -4,6 +4,7 @@
 
 from imod.mf6.chd import ConstantHead
 from imod.mf6.dis import StructuredDiscretization
+from imod.mf6.disu import LowLevelUnstructuredDiscretization
 from imod.mf6.disv import VerticesDiscretization
 from imod.mf6.drn import Drainage
 from imod.mf6.evt import Evapotranspiration
@@ -18,7 +19,7 @@
 from imod.mf6.model import GroundwaterFlowModel
 from imod.mf6.npf import NodePropertyFlow
 from imod.mf6.oc import OutputControl
-from imod.mf6.out import open_cbc, open_hds, read_cbc_headers, read_grb
+from imod.mf6.out import open_cbc, open_hds, open_hds_like, read_cbc_headers, read_grb
 from imod.mf6.rch import Recharge
 from imod.mf6.riv import River
 from imod.mf6.simulation import Modflow6Simulation
diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
index 68db74130..9e4a1912e 100644
--- a/imod/mf6/disu.py
+++ b/imod/mf6/disu.py
@@ -66,28 +66,29 @@ class LowLevelUnstructuredDiscretization(Package):
     ----------
     xorigin: float
     yorigin: float
-    top: array of floats (xr.DataArray)
-    bottom: array of floats (xr.DataArray)
-    area: array of floats (xr.DataArray)
-    iac: array of integers
-    ja: array of integers
-    ihc: array of integers
-    cl12: array of floats
-    hwva: array of floats
+    top: xr.DataArray of floats
+    bot: xr.DataArray of floats
+    area: xr.DataArray of floats
+    iac: xr.DataArray of integers
+    ja: xr.DataArray of integers
+    ihc: xr.DataArray of integers
+    cl12: xr.DataArray of floats
+    hwva: xr.DataArray of floats
     """
 
     _pkg_id = "disu"
     _grid_data = {
         "top": np.float64,
-        "bottom": np.float64,
+        "bot": np.float64,
         "area": np.float64,
         "iac": np.int32,
         "ja": np.int32,
         "ihc": np.int32,
         "cl12": np.float64,
         "hwva": np.float64,
+        "idomain": np.int32,
     }
-    _keyword_map = {"bottom": "bot"}
+    _keyword_map = {}
     _template = Package._initialize_template(_pkg_id)
 
     def __init__(
@@ -95,25 +96,28 @@ def __init__(
         xorigin,
         yorigin,
         top,
-        bottom,
+        bot,
         area,
         iac,
         ja,
         ihc,
         cl12,
         hwva,
+        idomain=None,
     ):
         super().__init__(locals())
         self.dataset["xorigin"] = xorigin
         self.dataset["yorigin"] = yorigin
         self.dataset["top"] = top
-        self.dataset["bottom"] = bottom
+        self.dataset["bot"] = bot
         self.dataset["area"] = area
         self.dataset["iac"] = iac
         self.dataset["ja"] = ja
         self.dataset["ihc"] = ihc
         self.dataset["cl12"] = cl12
         self.dataset["hwva"] = hwva
+        if idomain is not None:
+            self.dataset["idomain"] = idomain
 
     def render(self, directory, pkgname, globaltimes, binary):
         disdirectory = directory / pkgname
@@ -126,45 +130,38 @@ def render(self, directory, pkgname, globaltimes, binary):
         d["nja"] = int(self.dataset["iac"].sum())
 
         # Grid data
-        d["top"] = self._compose_values(
-            self.dataset["top"], disdirectory, "top", binary=binary
-        )[1][0]
-        d["bot"] = self._compose_values(
-            self["bottom"], disdirectory, "bot", binary=binary
-        )[1][0]
-        d["area"] = self._compose_values(
-            self["area"], disdirectory, "area", binary=binary
-        )[1][0]
-
-        # Connection data
-        d["iac"] = self._compose_values(
-            self["iac"], disdirectory, "iac", binary=binary
-        )[1][0]
-        d["ja"] = self._compose_values(self["ja"], disdirectory, "ja", binary=binary)[
-            1
-        ][0]
-        d["ihc"] = self._compose_values(
-            self["ihc"], disdirectory, "ihc", binary=binary
-        )[1][0]
-        d["cl12"] = self._compose_values(
-            self["cl12"], disdirectory, "cl12", binary=binary
-        )[1][0]
-        d["hwva"] = self._compose_values(
-            self["hwva"], disdirectory, "hwva", binary=binary
-        )[1][0]
+        for varname in self._grid_data:
+            if varname in self.dataset:
+                key = self._keyword_map.get(varname, varname)
+                d[varname] = self._compose_values(
+                    self.dataset[varname], disdirectory, key, binary=binary
+                )[1][0]
 
         return self._template.render(d)
 
     @staticmethod
-    def from_structured(
+    def from_dis(
         top,
         bottom,
         idomain,
+        reduce_nodes=False,
     ):
+        """
+        Parameters
+        ----------
+        reduce_nodes: bool, optional. Default: False.
+            Reduces the node numbering, discards cells when idomain <= 0.
+
+        Returns
+        -------
+        disu: LowLevelUnstructuredDiscretization
+        cell_ids: ndarray of integers.
+            Only provided if ``reduce_nodes`` is ``True``.
+        """
         x = idomain.coords["x"]
         y = idomain.coords["y"]
         layer = idomain.coords["layer"]
-        active = idomain.values > 0
+        active = idomain.values.ravel() > 0
 
         ncolumn = x.size
         nrow = y.size
@@ -174,17 +171,19 @@ def from_structured(
         dx, xmin, _ = imod.util.coord_reference(x)
         dy, ymin, _ = imod.util.coord_reference(y)
 
-        # MODFLOW6 expects the ja values to contain the cell number first
-        # while the row should be otherwise sorted ascending.
-        # scipy.spare.csr_matrix will sort the values ascending, but
-        # would not put the cell number first. To ensure this, we use
-        # the values as well as i and j; we sort on the zeros (thereby ensuring
-        # it results as a first value per column), but the actual value
-        # is the (negative) cell number (in v).
+        # MODFLOW6 expects the ja values to contain the cell number first while
+        # the row should be otherwise sorted ascending. scipy.spare.csr_matrix
+        # will sort the values ascending, but would not put the cell number
+        # first. To ensure this, we use the values as well as i and j; we sort
+        # on the zeros (thereby ensuring it results as a first value per
+        # column), but the actual value is the (negative) cell number (in v).
         ii, jj = _structured_connectivity(idomain)
         ii += 1
         jj += 1
-        nodes = np.arange(1, size + 1, dtype=np.int32)[active.ravel()]
+        nodes = np.arange(1, size + 1, dtype=np.int32)
+        if reduce_nodes:
+            nodes = nodes[active.ravel()]
+
         zeros = np.zeros_like(nodes)
         i = np.concatenate([nodes, ii, jj])
         j = np.concatenate([zeros, jj, ii])
@@ -194,16 +193,25 @@ def from_structured(
         # This entry does not require data in ihc, cl12, hwva.
         is_node = csr.data < 0
 
-        # Constructing the CSR matrix will have sorted all the values are
-        # required by MODFLOW6. However, we're using the original structured
-        # numbering, which includes inactive cells.
-        # For MODFLOW6, we use the reduced numbering, excluding all inactive
-        # cells. This means getting rid of empty rows (iac), generating (via
-        # cumsum) new numbers, and extracting them in the right order.
         nnz = csr.getnnz(axis=1)
-        iac = nnz[nnz > 0]
-        ja_index = np.abs(csr.data) - 1  # Get rid of negative values temporarily.
-        ja = active.ravel().cumsum()[ja_index]
+        if reduce_nodes:
+            # Constructing the CSR matrix will have sorted all the values are
+            # required by MODFLOW6. However, we're using the original structured
+            # numbering, which includes inactive cells.
+            # For MODFLOW6, we use the reduced numbering if reduce_nodes is True,
+            # excluding all inactive cells. This means getting rid of empty rows
+            # (iac), generating (via cumsum) new numbers, and extracting them in
+            # the right order.
+            iac = nnz[nnz > 0]
+            ja_index = np.abs(csr.data) - 1  # Get rid of negative values temporarily.
+            ja = active.cumsum()[ja_index]
+        else:
+            # In this case, inactive cells are included as well. They have no
+            # connections to other cells and form empty rows (0 in iac), but
+            # are still included. There is no need to update the cell numbers
+            # in this case.
+            iac = nnz[1:]  # Cell 0 does not exist.
+            ja = csr.data
 
         # From CSR back to COO form
         # connectivity for every cell: n -> m
@@ -253,15 +261,36 @@ def from_structured(
 
         # Set "node" and "nja" as the dimension in accordance with MODFLOW6.
         # Should probably be updated if we could move to UGRID3D...
-        return LowLevelUnstructuredDiscretization(
-            xorigin=xmin,
-            yorigin=ymin,
-            top=xr.DataArray(top.values[active], dims=["node"]),
-            bottom=xr.DataArray(bottom.values[active], dims=["node"]),
-            area=xr.DataArray(area[active.ravel()], dims=["node"]),
-            iac=xr.DataArray(iac, dims=["node"]),
-            ja=xr.DataArray(ja, dims=["nja"]),
-            ihc=xr.DataArray(ihc, dims=["nja"]),
-            cl12=xr.DataArray(cl12, dims=["nja"]),
-            hwva=xr.DataArray(hwva, dims=["nja"]),
-        )
+        if reduce_nodes:
+            # If we reduce nodes, we should only take active cells from top,
+            # bottom, area. There is no need to include an idomain: all defined
+            # cells are active.
+            disu = LowLevelUnstructuredDiscretization(
+                xorigin=xmin,
+                yorigin=ymin,
+                top=xr.DataArray(top.values.ravel()[active], dims=["node"]),
+                bot=xr.DataArray(bottom.values.ravel()[active], dims=["node"]),
+                area=xr.DataArray(area[active], dims=["node"]),
+                iac=xr.DataArray(iac, dims=["node"]),
+                ja=xr.DataArray(ja, dims=["nja"]),
+                ihc=xr.DataArray(ihc, dims=["nja"]),
+                cl12=xr.DataArray(cl12, dims=["nja"]),
+                hwva=xr.DataArray(hwva, dims=["nja"]),
+            )
+            cell_ids = np.cumsum(active) - 1
+            cell_ids[~active] = -1
+            return disu, cell_ids
+        else:
+            return LowLevelUnstructuredDiscretization(
+                xorigin=xmin,
+                yorigin=ymin,
+                top=xr.DataArray(top.values.ravel(), dims=["node"]),
+                bot=xr.DataArray(bottom.values.ravel(), dims=["node"]),
+                area=xr.DataArray(area, dims=["node"]),
+                iac=xr.DataArray(iac, dims=["node"]),
+                ja=xr.DataArray(ja, dims=["nja"]),
+                ihc=xr.DataArray(ihc, dims=["nja"]),
+                cl12=xr.DataArray(cl12, dims=["nja"]),
+                hwva=xr.DataArray(hwva, dims=["nja"]),
+                idomain=xr.DataArray(active.astype(np.int32), dims=["node"]),
+            )
diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index 0960066d4..c9da0feda 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -57,7 +57,8 @@ def disu_recarr(arrdict, layer, notnull):
     # Initialize the structured array
     nrow = notnull.sum()
     recarr = np.empty(nrow, dtype=sparse_dtype)
-    recarr["node"] = np.argwhere(notnull) + 1
+    # Argwhere returns an index_array with dims (N, a.ndims)
+    recarr["node"] = (np.argwhere(notnull) + 1)[:, 0]
     return recarr
 
 
@@ -369,14 +370,15 @@ def write(self, directory, pkgname, globaltimes, binary):
                 pkgdirectory.mkdir(exist_ok=True, parents=True)
                 for varname, dtype in self._grid_data.items():
                     key = self._keyword_map.get(varname, varname)
-                    da = self.dataset[varname]
-                    if self._is_xy_data(da):
-                        if binary:
-                            path = pkgdirectory / f"{key}.bin"
-                            self.write_binary_griddata(path, da, dtype)
-                        else:
-                            path = pkgdirectory / f"{key}.dat"
-                            self.write_text_griddata(path, da, dtype)
+                    if varname in self.dataset:
+                        da = self.dataset[varname]
+                        if self._is_xy_data(da):
+                            if binary:
+                                path = pkgdirectory / f"{key}.bin"
+                                self.write_binary_griddata(path, da, dtype)
+                            else:
+                                path = pkgdirectory / f"{key}.dat"
+                                self.write_text_griddata(path, da, dtype)
 
     @classmethod
     def filter_and_rename(cls, content: Dict[str, Any]) -> Dict[str, Any]:
@@ -472,17 +474,50 @@ def write_netcdf(self, directory, pkgname, aggregate_layers=False):
         spatial_ds.to_netcdf(path)
         return has_dims
 
-    def _to_disu(self, numbers):
-        structured = self.dataset.expand_dims("layer")
+    def _dis_to_disu(self, cell_ids=None):
+        structured = self.dataset
+        if "layer" not in structured.coords:
+            raise ValueError("layer coordinate is required")
+        if "layer" not in structured.dims:
+            structured = self.dataset.expand_dims("layer")
+
         # Stack will automatically broadcast to layer
         dataset = structured.stack(node=("layer", "y", "x"))
         layers = structured["layer"].values
         ncell_per_layer = structured["y"].size * structured["x"].size
         offset = (layers - 1) * ncell_per_layer
-        index = np.add.outer(offset, np.arange(ncell_per_layer))
-        dataset = dataset.assign_coords(node=numbers[index])
+        index = np.add.outer(offset, np.arange(ncell_per_layer)).ravel()
+
+        if cell_ids is not None:
+            node = cell_ids[index]
+        else:
+            node = index
+
+        dataset = dataset.assign_coords(node=node)
         return self.__class__(**dataset)
 
+    def to_disu(self, cell_ids=None):
+        """
+        Parameters
+        ----------
+        cell_ids: np.ndarray of integers, optional.
+            DISU cell IDs. Should contain values of -1 for inactive cells.
+
+        Returns
+        -------
+        disu_package: Any
+            Package in DISU form.
+        """
+        if isinstance(self.dataset, xr.Dataset):
+            return self._dis_to_disu(cell_ids)
+        elif isinstance(self.dataset.xu.UgridDataset):
+            raise NotImplementedError
+        else:
+            raise TypeError(
+                "Expected xarray.Dataset or xugrid.UgridDataset. "
+                f"Found: {type(self.dataset)}"
+            )
+
 
 class BoundaryCondition(Package, abc.ABC):
     """
diff --git a/imod/templates/mf6/gwf-disu.j2 b/imod/templates/mf6/gwf-disu.j2
index f1c2d7576..2990e1aa9 100644
--- a/imod/templates/mf6/gwf-disu.j2
+++ b/imod/templates/mf6/gwf-disu.j2
@@ -25,6 +25,9 @@ begin griddata
     {{bot}}
   area
     {{area}}
+{% if idomain is defined %}  idomain
+    {{idomain}}
+{% endif -%}
 end griddata
 
 begin connectiondata

From ae7c93dc5bdceec410daa336c964d066a256e650 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:54:21 +0200
Subject: [PATCH 10/35] newline in error message

---
 imod/mf6/read_input/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index 8e9df392c..b8f5374a2 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -181,7 +181,7 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
         warnings.warn(
             "Cannot convert time coordinate to datetime."
             "Falling back to (unitless) floating point time coordinate.\n"
-            f"time_units should be one of: {possibilities}; "
+            f"time_units should be one of: {possibilities}.\n"
             "start_date_time should be according to ISO 8601."
         )
         times = np.full(cumulative_length.size, 0.0)

From ed7a1621f56fa0954d6dce4de4e9a9c47424cee4 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:54:33 +0200
Subject: [PATCH 11/35] Add a data array reshape in utils

---
 imod/util.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/imod/util.py b/imod/util.py
index bf6d6cf04..2d069343b 100644
--- a/imod/util.py
+++ b/imod/util.py
@@ -1073,6 +1073,38 @@ def _replace(
     )
 
 
+def reshape(da: xr.DataArray, coords: Dict[str, Any]):
+    """
+    Gives a new shape and coordinates to a DataArray without changing its data.
+
+    Parameters
+    ----------
+    da: xr.DataArray
+    coords: dict
+        Coordinates of the reshaped DataArray. The key order determines
+        dimension order of the reshaped array. The sizes of the values
+        determines the size of the dimensions. Scalar values are expanded into
+        size 1 arrays. Multi-dimensional coordinate values are not allowed.
+
+    Returns
+    -------
+    reshaped: xr.DataArray
+
+    See also
+    --------
+    DataArray.stack
+    DataArray.unstack
+    """
+    coords = {k: np.atleast_1d(v) for k, v in coords.items()}
+    shape = [v.size for v in coords.values()]
+    for key, value in coords.items():
+        if value.ndim != 1:
+            raise ValueError(
+                f"{key} has ndim {value.ndim}. Can only reshape with 1D coordinates."
+            )
+    return xr.DataArray(da.data.reshape(shape), coords=coords, dims=list(coords.keys()))
+
+
 class MissingOptionalModule:
     """
     Presents a clear error for optional modules.

From f198ee734336058376a7881895bb6b6239247247 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:55:07 +0200
Subject: [PATCH 12/35] More comments on mf6 header sizes in output files

---
 imod/mf6/out/dis.py  | 9 +++++----
 imod/mf6/out/disv.py | 2 +-
 2 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/imod/mf6/out/dis.py b/imod/mf6/out/dis.py
index 1a6dea3ef..6ea5f7a9f 100644
--- a/imod/mf6/out/dis.py
+++ b/imod/mf6/out/dis.py
@@ -74,11 +74,11 @@ def read_times(
     """
     times = np.empty(ntime, dtype=np.float64)
 
-    # Compute how much to skip to the next timestamp
-    start_of_header = 16
-    rest_of_header = 28
+    # Compute how much to skip to the next timestamp.
+    start_of_header = 16  # KSTP(4), KPER(4), PERTIM(8)
+    rest_of_header = 28  # TEXT(16), NCOL(4), NROW(4), ILAY(4)
     data_single_layer = nrow * ncol * 8
-    header = 52
+    header = 52  # start_of_header + TOTIM(8) + rest_of_header
     nskip = (
         rest_of_header
         + data_single_layer
@@ -117,6 +117,7 @@ def read_hds_timestep(
 def open_hds(path: FilePath, d: Dict[str, Any], dry_nan: bool) -> xr.DataArray:
     nlayer, nrow, ncol = d["nlayer"], d["nrow"], d["ncol"]
     filesize = os.path.getsize(path)
+    # 52 is header size; 8 is size of double.
     ntime = filesize // (nlayer * (52 + (nrow * ncol * 8)))
     times = read_times(path, ntime, nlayer, nrow, ncol)
     coords = d["coords"]
diff --git a/imod/mf6/out/disv.py b/imod/mf6/out/disv.py
index d53109b98..5093a348d 100644
--- a/imod/mf6/out/disv.py
+++ b/imod/mf6/out/disv.py
@@ -130,7 +130,7 @@ def read_hds_timestep(
         f.seek(pos)
         a1d = np.empty(nlayer * ncells_per_layer, dtype=np.float64)
         for k in range(nlayer):
-            f.seek(52, 1)  # skip kstp, kper, pertime
+            f.seek(52, 1)  # skip header
             a1d[k * ncells_per_layer : (k + 1) * ncells_per_layer] = np.fromfile(
                 f, np.float64, ncells_per_layer
             )

From e6d54a3a184d63887c6cb87b51e0048f22c4abb2 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:55:20 +0200
Subject: [PATCH 13/35] Read DISU head output

---
 imod/mf6/out/disu.py | 51 ++++++++++++++++++++++++++++++++++++--------
 1 file changed, 42 insertions(+), 9 deletions(-)

diff --git a/imod/mf6/out/disu.py b/imod/mf6/out/disu.py
index f7c0372d9..f884d8157 100644
--- a/imod/mf6/out/disu.py
+++ b/imod/mf6/out/disu.py
@@ -1,12 +1,14 @@
+import os
 import struct
 from typing import Any, BinaryIO, Dict, List
 
+import dask
 import numpy as np
 import xarray as xr
 import xugrid as xu
 
 from . import cbc
-from .common import FilePath, _grb_text
+from .common import FilePath, _grb_text, _to_nan
 
 
 def read_grb(f: BinaryIO, ntxt: int, lentxt: int) -> Dict[str, Any]:
@@ -58,17 +60,48 @@ def read_grb(f: BinaryIO, ntxt: int, lentxt: int) -> Dict[str, Any]:
 
 
 def read_times(path: FilePath, ntime: int, ncell: int):
-    raise NotImplementedError
-
-
-def read_hds_timestep(
-    path: FilePath, nlayer: int, ncell_per_layer: int, dry_nan: bool, pos: int
-):
-    raise NotImplementedError
+    """
+    Reads all total simulation times.
+    """
+    times = np.empty(ntime, dtype=np.float64)
+
+    # Compute how much to skip to the next timestamp
+    start_of_header = 16  # KSTP(4), KPER(4), PERTIM(8)
+    rest_of_header = 28  # TEXT(16), NCOL(4), NROW(4), ILAY(4)
+    data = ncell * 8
+    nskip = rest_of_header + data + start_of_header
+    with open(path, "rb") as f:
+        f.seek(start_of_header)
+        for i in range(ntime):
+            times[i] = struct.unpack("d", f.read(8))[0]
+            f.seek(nskip, 1)
+    return times
+
+
+def read_hds_timestep(path: FilePath, ncell: int, dry_nan: bool, pos: int):
+    with open(path, "rb") as f:
+        f.seek(pos + 52)
+        a1d = np.fromfile(f, np.float64, ncell)
+    return _to_nan(a1d, dry_nan)
 
 
 def open_hds(path: FilePath, d: Dict[str, Any], dry_nan: bool) -> xr.DataArray:
-    raise NotImplementedError
+    ncell = d["ncell"]
+    filesize = os.path.getsize(path)
+    ntime = filesize // (52 + ncell * 8)
+    times = read_times(path, ntime, ncell)
+    coords = d["coords"]
+    coords["time"] = times
+
+    dask_list = []
+    for i in range(ntime):
+        pos = i * (52 + ncell * 8)
+        a = dask.delayed(read_hds_timestep)(path, ncell, dry_nan, pos)
+        x = dask.array.from_delayed(a, shape=(ncell,), dtype=np.float64)
+        dask_list.append(x)
+
+    daskarr = dask.array.stack(dask_list, axis=0)
+    return xr.DataArray(daskarr, coords, ("time", "node"), name=d["name"])
 
 
 def open_imeth1_budgets(

From 724cbde4ffde4688c6c4d985dd45e9c67b9c6610 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 10 May 2022 10:55:38 +0200
Subject: [PATCH 14/35] Add name to open_hds_like

---
 imod/mf6/out/__init__.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/imod/mf6/out/__init__.py b/imod/mf6/out/__init__.py
index 1ba37bc5f..3807e58d7 100644
--- a/imod/mf6/out/__init__.py
+++ b/imod/mf6/out/__init__.py
@@ -151,10 +151,12 @@ def open_hds_like(
     # TODO: check shape with hds metadata.
     if isinstance(like, xr.DataArray):
         d = dis.grid_info(like)
+        d["name"] = "head"
         return dis.open_hds(path, d, dry_nan)
 
     elif isinstance(like, xu.UgridDataArray):
         d = disv.grid_info(like)
+        d["name"] = "head"
         return disv.open_hds(path, d, dry_nan)
 
     else:

From 2b4819d69caf6ef2aade6facfb2d50af1ab01c5e Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 16 May 2022 11:06:14 +0200
Subject: [PATCH 15/35] Uncomment jit statement

---
 imod/mf6/disu.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
index 9e4a1912e..a58ecfdda 100644
--- a/imod/mf6/disu.py
+++ b/imod/mf6/disu.py
@@ -14,8 +14,8 @@ def _number(k, i, j, nrow, ncolumn):
     return k * (nrow * ncolumn) + i * ncolumn + j
 
 
-# @nb.njit
-def _structured_connectivity(idomain):
+@nb.njit
+def _structured_connectivity(idomain: IntArray):
     nlayer, nrow, ncolumn = idomain.shape
     # Pre-allocate: structured connectivity implies maximum of 8 neighbors
     nconnection = idomain.size * 8
@@ -177,7 +177,7 @@ def from_dis(
         # first. To ensure this, we use the values as well as i and j; we sort
         # on the zeros (thereby ensuring it results as a first value per
         # column), but the actual value is the (negative) cell number (in v).
-        ii, jj = _structured_connectivity(idomain)
+        ii, jj = _structured_connectivity(idomain.values)
         ii += 1
         jj += 1
         nodes = np.arange(1, size + 1, dtype=np.int32)

From 25cb83d68bd5e4ec18ddc1ff10609c85d676772f Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 16 May 2022 11:06:29 +0200
Subject: [PATCH 16/35] Warning message formatting

---
 imod/mf6/read_input/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index b8f5374a2..e87b75d5b 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -180,7 +180,7 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
         possibilities = ", ".join(list(TIME_UNITS.keys())[1:])
         warnings.warn(
             "Cannot convert time coordinate to datetime."
-            "Falling back to (unitless) floating point time coordinate.\n"
+            "Falling back to (unitless) floating point time coordinate. \n"
             f"time_units should be one of: {possibilities}.\n"
             "start_date_time should be according to ISO 8601."
         )

From fb504bcf9942e1be231624cf31f503042934e4cc Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 18 May 2022 18:05:35 +0200
Subject: [PATCH 17/35] Add tests for read_input/common.py and grid_data.py

---
 imod/mf6/read_input/common.py                 |  35 +--
 imod/mf6/read_input/grid_data.py              |   2 +-
 .../test_mf6_input_reading_common.py          | 247 ++++++++++++++++++
 .../test_mf6_input_reading_griddata.py        | 121 +++++++++
 4 files changed, 387 insertions(+), 18 deletions(-)
 create mode 100644 tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
 create mode 100644 tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py

diff --git a/imod/mf6/read_input/common.py b/imod/mf6/read_input/common.py
index 1142b5149..5bb50c6c1 100644
--- a/imod/mf6/read_input/common.py
+++ b/imod/mf6/read_input/common.py
@@ -27,7 +27,8 @@ def flatten(iterable: Any) -> List[Any]:
 def split_line(line: str) -> List[str]:
     # Maybe check: https://stackoverflow.com/questions/36165050/python-equivalent-of-fortran-list-directed-input
     # Split on comma and whitespace, like a FORTRAN read would do.
-    return flatten([part.split(",") for part in line.split()])
+    flat = flatten([part.split(",") for part in line.split()])
+    return [part for part in flat if part != ""]
 
 
 def to_float(string: str) -> float:
@@ -137,6 +138,22 @@ def advance_to_header(f: IO[str], section) -> None:
     raise ValueError(f'"{start}" is not present in file {f.name}')
 
 
+def parse_option(stripped: str, fname: str) -> Tuple[str, Any]:
+    separated = stripped.split()
+    nwords = len(separated)
+    if nwords == 0:
+        raise ValueError(f"Cannot parse option in {fname}")
+    elif nwords == 1:
+        key = separated[0]
+        value = True
+    elif nwords == 2:
+        key, value = separated
+    else:
+        key = separated[0]
+        value = separated[1:]
+    return key, value
+
+
 def read_key_value_block(f: IO[str], parse: Callable) -> Dict[str, str]:
     fname = f.name
     content = {}
@@ -180,22 +197,6 @@ def read_iterable_block(f: IO[str], parse: Callable) -> List[Any]:
     raise ValueError(f'"end" of block is not present in file {fname}')
 
 
-def parse_option(stripped: str, fname: str) -> Tuple[str, Any]:
-    separated = stripped.split()
-    nwords = len(separated)
-    if nwords == 0:
-        raise ValueError(f"Cannot parse {stripped} in {fname}")
-    elif nwords == 1:
-        key = separated[0]
-        value = True
-    elif nwords == 2:
-        key, value = separated
-    else:
-        key = separated[0]
-        value = separated[1:]
-    return key, value
-
-
 def parse_dimension(stripped: str, fname: str) -> Tuple[str, int]:
     key, value = parse_option(stripped, fname)
     return key, int(value)
diff --git a/imod/mf6/read_input/grid_data.py b/imod/mf6/read_input/grid_data.py
index 4a0780b70..bee18fdba 100644
--- a/imod/mf6/read_input/grid_data.py
+++ b/imod/mf6/read_input/grid_data.py
@@ -167,7 +167,7 @@ def read_array(
         array = dask.array.from_delayed(a, shape=shape, dtype=dtype)
     else:
         raise ValueError(
-            'Expected "constant", "internal" or "open/close". Received instead: '
+            'Expected "constant", "internal", or "open/close". Received instead: '
             f"{stripped}"
         )
 
diff --git a/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py b/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
new file mode 100644
index 000000000..b2b9328b9
--- /dev/null
+++ b/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
@@ -0,0 +1,247 @@
+import numpy as np
+import pytest
+
+from imod.mf6.read_input import common as cm
+
+
+def isa(value, kind, expected):
+    return isinstance(value, kind) and value == expected
+
+
+def test_end_of_file():
+    assert cm.end_of_file("") is True
+    assert cm.end_of_file(" ") is False
+    assert cm.end_of_file("\n") is False
+
+
+def test_strip_line():
+    assert cm.strip_line("abc") == "abc"
+    assert cm.strip_line("abc ") == "abc"
+    assert cm.strip_line("ABC ") == "abc"
+    assert cm.strip_line("abc#def") == "abc"
+    assert cm.strip_line("abc #def") == "abc"
+    assert cm.strip_line(" abc ##def") == "abc"
+
+
+def test_flatten():
+    assert cm.flatten([[1, 2], [3, 4]]) == [1, 2, 3, 4]
+    assert cm.flatten([["abc", "def"], ["ghi"]]) == ["abc", "def", "ghi"]
+
+
+def test_split_line():
+    assert cm.split_line("abc def ghi") == ["abc", "def", "ghi"]
+    assert cm.split_line("abc,def,ghi") == ["abc", "def", "ghi"]
+    assert cm.split_line("abc,def, ghi") == ["abc", "def", "ghi"]
+    assert cm.split_line("abc ,def, ghi ") == ["abc", "def", "ghi"]
+
+
+def test_to_float():
+    assert isa(cm.to_float("1"), float, 1.0)
+    assert isa(cm.to_float("1.0"), float, 1.0)
+    assert isa(cm.to_float("1.0d0"), float, 1.0)
+    assert isa(cm.to_float("1.0e0"), float, 1.0)
+    assert isa(cm.to_float("1.0+0"), float, 1.0)
+    assert isa(cm.to_float("1.0-0"), float, 1.0)
+    assert isa(cm.to_float("1.0e-0"), float, 1.0)
+    assert isa(cm.to_float("1.0e+0"), float, 1.0)
+
+
+def test_find_entry():
+    assert isa(cm.find_entry("abc factor 1", "factor", float), float, 1.0)
+    assert isa(cm.find_entry("abc factor 1", "factor", int), int, 1)
+    assert isa(cm.find_entry("abc factor 1", "factor", str), str, "1")
+    assert cm.find_entry("abc 1", "factor", float) is None
+
+
+def test_read_internal(tmp_path):
+    path1 = tmp_path / "internal-1.dat"
+    with open(path1, "w") as f:
+        f.write("1 2 3 4")
+    # max_rows should read all lines, unless max_rows is exceeded.
+    with open(path1) as f:
+        a = cm.read_internal(f, int, 1)
+    with open(path1) as f:
+        b = cm.read_internal(f, int, 2)
+    assert np.issubdtype(a.dtype, np.integer)
+    assert np.array_equal(a, b)
+
+    path2 = tmp_path / "internal-2.dat"
+    with open(path2, "w") as f:
+        f.write("1 2\n3 4")
+    with open(path2) as f:
+        a = cm.read_internal(f, float, 1)
+    with open(path2) as f:
+        b = cm.read_internal(f, int, 2)
+    assert np.issubdtype(a.dtype, np.floating)
+    assert np.issubdtype(b.dtype, np.integer)
+    assert a.size == 2
+    assert b.size == 4
+
+
+def test_read_external_binaryfile(tmp_path):
+    path1 = tmp_path / "external-1.bin"
+    a = np.ones((5, 5))
+    a.tofile(path1)
+    b = cm.read_external_binaryfile(path1, np.float64, 25)
+    assert b.shape == (25,)
+    b = cm.read_external_binaryfile(path1, np.float64, 10)
+    assert b.shape == (10,)
+
+    dtype = np.dtype([("node", np.int32), ("value", np.float32)])
+    a = np.ones((5,), dtype)
+    path2 = tmp_path / "external-2.bin"
+    a.tofile(path2)
+    b = cm.read_external_binaryfile(path2, dtype, 5)
+    assert np.array_equal(a, b)
+
+
+def test_read_fortran_deflated_text_array(tmp_path):
+    path1 = tmp_path / "deflated-1.txt"
+    with open(path1, "w") as f:
+        f.write("1.0\n2*2.0\n3*3.0")
+    a = cm.read_fortran_deflated_text_array(path1, float, 3)
+    b = np.array([1.0, 2.0, 2.0])
+    assert np.array_equal(a, b)
+
+    a = cm.read_fortran_deflated_text_array(path1, float, 6)
+    b = np.array([1.0, 2.0, 2.0, 3.0, 3.0, 3.0])
+    assert np.array_equal(a, b)
+
+
+def test_read_external_textfile(tmp_path):
+    path1 = tmp_path / "external.dat"
+    with open(path1, "w") as f:
+        f.write("1.0 2.0 3.0")
+    a = cm.read_external_textfile(path1, float, 6)
+    b = np.array([1.0, 2.0, 3.0])
+    assert np.array_equal(a, b)
+
+    path2 = tmp_path / "deflated-1.txt"
+    with open(path2, "w") as f:
+        f.write("1.0\n2*2.0\n3*3.0")
+    a = cm.read_external_textfile(path2, float, 6)
+    b = np.array([1.0, 2.0, 2.0, 3.0, 3.0, 3.0])
+    assert np.array_equal(a, b)
+
+
+def test_advance_to_header(tmp_path):
+    path = tmp_path / "header.txt"
+    content = "\n".join(
+        [
+            "",
+            "begin options",
+            "1",
+            "end options\n",
+            "\n",
+            "begin griddata",
+            "2",
+            "end griddata",
+        ]
+    )
+    with open(path, "w") as f:
+        f.write(content)
+
+    with open(path) as f:
+        cm.advance_to_header(f, "options")
+        assert f.readline() == "1\n"
+        f.readline()  # read "end options" line
+        cm.advance_to_header(f, "griddata")
+        assert f.readline() == "2\n"
+
+        with pytest.raises(ValueError, match='"begin perioddata" is not present'):
+            cm.advance_to_header(f, "perioddata")
+
+
+def test_parse_option():
+    with pytest.raises(ValueError, match="Cannot parse option in options.txt"):
+        cm.parse_option("", "options.txt")
+    assert cm.parse_option("save_flows", "options.txt") == ("save_flows", True)
+    assert cm.parse_option("variablecv dewatered", "options.txt") == (
+        "variablecv",
+        "dewatered",
+    )
+    assert cm.parse_option("multiple a b c", "options.txt") == (
+        "multiple",
+        ["a", "b", "c"],
+    )
+
+
+def test_read_key_value_block(tmp_path):
+    path = tmp_path / "values.txt"
+    content = "\n".join(["", "save_flows", "variablecv dewatered", "", "end options"])
+
+    with open(path, "w") as f:
+        f.write(content)
+
+    with open(path) as f:
+        d = cm.read_key_value_block(f, cm.parse_option)
+    assert d == {"save_flows": True, "variablecv": "dewatered"}
+
+    path2 = tmp_path / "values-unterminated.txt"
+    content = "\n".join(
+        [
+            "",
+            "save_flows",
+            "variablecv dewatered",
+            "",
+        ]
+    )
+
+    with open(path2, "w") as f:
+        f.write(content)
+
+    with open(path2) as f:
+        with pytest.raises(ValueError, match='"end" of block is not present in file'):
+            cm.read_key_value_block(f, cm.parse_option)
+
+
+def test_read_iterable_block(tmp_path):
+    def parse(line, _):  # second arg is normally file name
+        return line.split()
+
+    path = tmp_path / "iterable-values.txt"
+    content = "\n".join(["1.0 1 1.0", "1.0 1 1.0", "1.0 1 1.0", "end perioddata"])
+    with open(path, "w") as f:
+        f.write(content)
+
+    with open(path) as f:
+        back = cm.read_iterable_block(f, parse)
+    assert back == [
+        ["1.0", "1", "1.0"],
+        ["1.0", "1", "1.0"],
+        ["1.0", "1", "1.0"],
+    ]
+
+    path2 = tmp_path / "iterable-values-unterminated.txt"
+    content = "\n".join(
+        [
+            "1.0 1 1.0",
+            "1.0 1 1.0",
+            "1.0 1 1.0",
+        ]
+    )
+    with open(path2, "w") as f:
+        f.write(content)
+    with open(path2) as f:
+        with pytest.raises(ValueError, match='"end" of block is not present in file'):
+            back = cm.read_iterable_block(f, parse)
+
+
+def test_parse_dimension():
+    assert cm.parse_dimension("  nper 3  ", "fname") == ("nper", 3)
+
+
+def test_advance_to_period(tmp_path):
+    path = tmp_path / "period.txt"
+    content = "\n".join(
+        [
+            "" "  begin period 1",
+            "2",
+        ]
+    )
+    with open(path, "w") as f:
+        f.write(content)
+
+    with open(path) as f:
+        cm.advance_to_period(f)
+        assert f.readline() == "2"
diff --git a/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py b/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
new file mode 100644
index 000000000..d1961ba1b
--- /dev/null
+++ b/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
@@ -0,0 +1,121 @@
+import dask.array
+import numpy as np
+import pytest
+
+from imod.mf6.read_input import grid_data as gd
+
+
+def test_advance_to_griddata_section(tmp_path):
+    path = tmp_path / "griddata.txt"
+    content = "\n".join(
+        [
+            "",
+            " idomain",
+            "   open/close GWF_1/dis/idomain.dat",
+            " botm layered",
+            "",
+        ]
+    )
+
+    with open(path, "w") as f:
+        f.write(content)
+
+    with open(path) as f:
+        assert gd.advance_to_griddata_section(f) == ("idomain", False)
+        f.readline()
+        assert gd.advance_to_griddata_section(f) == ("botm", True)
+        with pytest.raises(ValueError, match="No end of griddata specified"):
+            gd.advance_to_griddata_section(f)
+
+
+def test_shape_to_max_rows():
+    assert gd.shape_to_max_rows(shape=(4, 3, 2), layered=False) == (12, (4, 3, 2))
+    assert gd.shape_to_max_rows(shape=(4, 3, 2), layered=True) == (3, (3, 2))
+    assert gd.shape_to_max_rows(shape=(3, 2), layered=False) == (3, (3, 2))
+    assert gd.shape_to_max_rows(shape=(3, 2), layered=True) == (1, (2,))
+    assert gd.shape_to_max_rows(shape=(6,), layered=False) == (1, (6,))
+    with pytest.raises(
+        ValueError, match="LAYERED section detected. DISU does not support LAYERED"
+    ):
+        gd.shape_to_max_rows(shape=(6,), layered=True)
+    with pytest.raises(ValueError, match="length of shape should be 1, 2, or 3"):
+        gd.shape_to_max_rows(shape=(5, 4, 3, 2), layered=True)
+
+
+def test_constant():
+    shape = (3, 4, 5)
+    a = gd.constant(2.0, shape, np.float64)
+    assert isinstance(a, dask.array.Array)
+    b = a.compute()
+    assert np.allclose(b, 2.0)
+
+
+def test_read_internal_griddata(tmp_path):
+    path = tmp_path / "internal.txt"
+    with open(path, "w") as f:
+        f.write("1 2 3 4\n5 6 7 8")
+
+    with open(path) as f:
+        a = gd.read_internal_griddata(f, np.int32, (2, 4), 2)
+
+    assert a.shape == (2, 4)
+    assert np.array_equal(a.ravel(), np.arange(1, 9))
+
+
+def test_read_external_griddata(tmp_path):
+    path1 = tmp_path / "external.dat"
+    path2 = tmp_path / "external.bin"
+
+    a = np.ones((4, 2))
+    a.tofile(path1, sep=" ")
+    a.tofile(path2)  # binary
+
+    b = gd.read_external_griddata(path1, np.float64, (4, 2), False)
+    c = gd.read_external_griddata(path2, np.float64, (4, 2), True)
+    assert np.array_equal(a, b)
+    assert np.array_equal(a, c)
+
+
+def test_read_array(tmp_path):
+    blockfile_path = tmp_path / "blockfile.txt"
+    external_path = tmp_path / "external.dat"
+    external_binary_path = tmp_path / "external.bin"
+
+    content = "\n".join(
+        [
+            "top",
+            "  constant 200.0",
+            "idomain",
+            "  open/close external.dat",
+            "botm",
+            "  open/close external.bin (binary)",
+            "abc",
+        ]
+    )
+    with open(blockfile_path, "w") as f:
+        f.write(content)
+
+    shape = (3, 4, 5)
+    a = np.ones(shape, dtype=np.float64)
+    a.tofile(external_path, sep=" ")
+    a.tofile(external_binary_path)
+
+    with open(blockfile_path) as f:
+        f.readline()
+        top = gd.read_array(f, tmp_path, np.float64, max_rows=12, shape=shape)
+        f.readline()
+        idomain = gd.read_array(f, tmp_path, np.float64, max_rows=12, shape=shape)
+        f.readline()
+        botm = gd.read_array(f, tmp_path, np.float64, max_rows=12, shape=shape)
+        with pytest.raises(
+            ValueError, match='Expected "constant", "internal", or "open/close"'
+        ):
+            gd.read_array(f, tmp_path, np.float64, max_rows=12, shape=shape)
+
+    for a in (top, idomain, botm):
+        assert isinstance(a, dask.array.Array)
+        assert a.shape == shape
+
+    assert np.allclose(top, 200.0)
+    assert np.allclose(idomain, 1.0)
+    assert np.allclose(botm, 1.0)

From b122897e3a9fcb6e72d205f86ecef4773b66340e Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 19 May 2022 16:47:58 +0200
Subject: [PATCH 18/35] format

---
 imod/mf6/pkgbase.py    | 5 ++++-
 imod/mf6/simulation.py | 2 +-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index 8c620f2b7..58559e780 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -1,8 +1,8 @@
 import abc
 import inspect
 import pathlib
-from typing import Any, Dict
 from dataclasses import dataclass
+from typing import Any, Dict
 
 import jinja2
 import numpy as np
@@ -61,6 +61,8 @@ def disu_recarr(arrdict, layer, notnull):
     # Argwhere returns an index_array with dims (N, a.ndims)
     recarr["node"] = (np.argwhere(notnull) + 1)[:, 0]
     return recarr
+
+
 @dataclass
 class VariableMetaData:
     """
@@ -529,6 +531,7 @@ def to_disu(self, cell_ids=None):
                 "Expected xarray.Dataset or xugrid.UgridDataset. "
                 f"Found: {type(self.dataset)}"
             )
+
     def _check_types(self):
         for varname, metadata in self._metadata_dict.items():
             expected_dtype = metadata.dtype
diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py
index c7a00ec29..b0be98591 100644
--- a/imod/mf6/simulation.py
+++ b/imod/mf6/simulation.py
@@ -1,8 +1,8 @@
 import collections
 import pathlib
 import subprocess
-from typing import Union
 import warnings
+from typing import Union
 
 import jinja2
 import numpy as np

From 92d10b27747795b69bac9ac7100b08dd08ec65fd Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 19 May 2022 16:48:42 +0200
Subject: [PATCH 19/35] Tests for mf6 read_input: grid_data and list_input

---
 imod/mf6/read_input/list_input.py             |   6 +-
 .../test_mf6_input_reading_griddata.py        |  54 ++++-
 .../test_mf6_input_reading_listinput.py       | 217 ++++++++++++++++++
 3 files changed, 275 insertions(+), 2 deletions(-)
 create mode 100644 imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py

diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
index e75a97b07..1307e6eb1 100644
--- a/imod/mf6/read_input/list_input.py
+++ b/imod/mf6/read_input/list_input.py
@@ -38,6 +38,7 @@ def read_text_listinput(
 ) -> np.ndarray:
     # This function is over three times faster than:
     # recarr = np.loadtxt(path, dtype, max_rows=max_rows)
+    # I guess MODFLOW6 will also accept comma delimited?
     d = {key: value[0] for key, value in dtype.fields.items()}
     df = pd.read_csv(
         path,
@@ -138,8 +139,11 @@ def read_listinput(
         x = dask.delayed(read_external_listinput, nout=nout)(
             path, dtype, index_columns, fields, shape, binary, max_rows
         )
+
+        fieldtypes = [dtype.fields[field][0] for field in fields]
         variable_values = [
-            dask.array.from_delayed(a, shape=shape, dtype=dtype) for a in x
+            dask.array.from_delayed(a, shape=shape, dtype=dt)
+            for a, dt in zip(x, fieldtypes)
         ]
     else:
         # Move file position back one line, and try reading internal values.
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
index d1961ba1b..f16e668b2 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
@@ -117,5 +117,57 @@ def test_read_array(tmp_path):
         assert a.shape == shape
 
     assert np.allclose(top, 200.0)
-    assert np.allclose(idomain, 1.0)
+    assert np.allclose(idomain, 1)
     assert np.allclose(botm, 1.0)
+
+
+def test_read_griddata(tmp_path):
+    blockfile_path = tmp_path / "blockfile.txt"
+    external_path = tmp_path / "external.dat"
+    external_binary_path = tmp_path / "external.bin"
+
+    content = "\n".join(
+        [
+            "top",
+            "  constant 200.0",
+            "idomain",
+            "  open/close external.dat",
+            "botm",
+            "  open/close external.bin (binary)",
+            "end",
+        ]
+    )
+    with open(blockfile_path, "w") as f:
+        f.write(content)
+
+    sections = {
+        "top": (np.float64, gd.shape_to_max_rows),
+        "idomain": (np.int32, gd.shape_to_max_rows),
+        "botm": (np.float64, gd.shape_to_max_rows),
+    }
+    shape = (3, 4, 5)
+    a = np.ones(shape, dtype=np.int32)
+    a.tofile(external_path, sep=" ")
+    b = np.ones(shape, dtype=np.float64)
+    b.tofile(external_binary_path)
+
+    with open(blockfile_path) as f:
+        d = gd.read_griddata(f, tmp_path, sections, shape)
+
+    for key in ("top", "idomain", "botm"):
+        assert key in d
+        a = d[key]
+        assert isinstance(a, dask.array.Array)
+        assert a.shape == shape
+
+    assert np.allclose(d["top"], 200.0)
+    assert np.allclose(d["idomain"], 1)
+    assert np.allclose(d["botm"], 1.0)
+
+    dummy_path = tmp_path / "dummy.txt"
+    with open(dummy_path, "w") as f:
+        f.write("\n")
+
+    with open(dummy_path) as f:
+        with pytest.raises(ValueError, match="Error reading"):
+            d = gd.read_griddata(f, tmp_path, sections, shape)
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
new file mode 100644
index 000000000..473430dcd
--- /dev/null
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
@@ -0,0 +1,217 @@
+import dask.array
+import numpy as np
+
+from imod.mf6.read_input import list_input as li
+
+DIS_DTYPE = np.dtype(
+    [
+        ("layer", np.int32),
+        ("row", np.int32),
+        ("column", np.int32),
+        ("head", np.float64),
+        ("conductance", np.float64),
+    ]
+)
+
+
+def test_recarr_to_dense__dis():
+    dtype = DIS_DTYPE
+    recarr = np.array(
+        [
+            (1, 1, 1, 1.0, 10.0),
+            (2, 1, 1, 2.0, 10.0),
+            (3, 1, 1, 3.0, 10.0),
+        ],
+        dtype=dtype,
+    )
+
+    variables = li.recarr_to_dense(
+        recarr,
+        index_columns=["layer", "row", "column"],
+        fields=["head", "conductance"],
+        shape=(3, 4, 5),
+    )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+    a, b = variables
+    assert np.isfinite(a).sum() == 3
+    assert np.isfinite(b).sum() == 3
+    assert a[0, 0, 0] == 1.0
+    assert a[1, 0, 0] == 2.0
+    assert a[2, 0, 0] == 3.0
+
+
+def test_recarr_to_dense__disv():
+    dtype = np.dtype(
+        [
+            ("layer", np.int32),
+            ("cell2d", np.int32),
+            ("head", np.float64),
+            ("conductance", np.float64),
+        ]
+    )
+    recarr = np.array(
+        [
+            (1, 1, 1.0, 10.0),
+            (2, 1, 2.0, 10.0),
+            (3, 1, 3.0, 10.0),
+        ],
+        dtype=dtype,
+    )
+    variables = li.recarr_to_dense(
+        recarr,
+        index_columns=["layer", "cell2d"],
+        fields=["head", "conductance"],
+        shape=(3, 20),
+    )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+    a, b = variables
+    assert np.isfinite(a).sum() == 3
+    assert np.isfinite(b).sum() == 3
+    assert a[0, 0] == 1.0
+    assert a[1, 0] == 2.0
+    assert a[2, 0] == 3.0
+
+
+def test_recarr_to_dense__disu():
+    dtype = np.dtype(
+        [
+            ("node", np.int32),
+            ("head", np.float64),
+            ("conductance", np.float64),
+        ]
+    )
+    recarr = np.array(
+        [
+            (1, 1.0, 10.0),
+            (21, 2.0, 10.0),
+            (41, 3.0, 10.0),
+        ],
+        dtype=dtype,
+    )
+    variables = li.recarr_to_dense(
+        recarr,
+        index_columns=["node"],
+        fields=["head", "conductance"],
+        shape=(60,),
+    )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+    a, b = variables
+    assert np.isfinite(a).sum() == 3
+    assert np.isfinite(b).sum() == 3
+    assert a[0] == 1.0
+    assert a[20] == 2.0
+    assert a[40] == 3.0
+
+
+def test_read_text_listinput(tmp_path):
+    dtype = DIS_DTYPE
+    path = tmp_path / "listinput.dat"
+    content = "\n".join(
+        [
+            "# layer row column head conductance",
+            "1 1 1 1.0 10.0",
+            "2 1 1 2.0 10.0",
+            "3 1 1 3.0 10.0",
+        ]
+    )
+
+    with open(path, "w") as f:
+        f.write(content)
+
+    # Test for internal input as well, for an already opened file.
+    with open(path) as f:
+        variables = li.read_internal_listinput(
+            f,
+            dtype,
+            index_columns=["layer", "row", "column"],
+            fields=["head", "conductance"],
+            max_rows=3,
+            shape=(3, 4, 5),
+        )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+
+    variables = li.read_external_listinput(
+        path,
+        dtype,
+        index_columns=["layer", "row", "column"],
+        fields=["head", "conductance"],
+        shape=(3, 4, 5),
+        binary=False,
+        max_rows=3,
+    )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+
+
+def test_read_binary_listinput(tmp_path):
+    path = tmp_path / "listinput.bin"
+    dtype = DIS_DTYPE
+    recarr = np.array(
+        [
+            (1, 1, 1, 1.0, 10.0),
+            (2, 1, 1, 2.0, 10.0),
+            (3, 1, 1, 3.0, 10.0),
+        ],
+        dtype=dtype,
+    )
+    recarr.tofile(path)
+
+    variables = li.read_external_listinput(
+        path,
+        dtype,
+        index_columns=["layer", "row", "column"],
+        fields=["head", "conductance"],
+        shape=(3, 4, 5),
+        binary=True,
+        max_rows=3,
+    )
+    assert isinstance(variables, list)
+    assert len(variables) == 2
+
+
+def test_read_listinput(tmp_path):
+    dtype = DIS_DTYPE
+    path = "package-binary.txt"
+    content = "\n".join(
+        [
+            "open/close listinput.bin (binary)",
+        ]
+    )
+    with open(path, "w") as f:
+        f.write(content)
+
+    binpath = tmp_path / "listinput.bin"
+    dtype = DIS_DTYPE
+    recarr = np.array(
+        [
+            (1, 1, 1, 1.0, 10.0),
+            (2, 1, 1, 2.0, 10.0),
+            (3, 1, 1, 3.0, 10.0),
+        ],
+        dtype=dtype,
+    )
+    recarr.tofile(binpath)
+
+    with open(path) as f:
+        variables = li.read_listinput(
+            f,
+            tmp_path,
+            dtype,
+            index_columns=["layer", "row", "column"],
+            fields=["head", "conductance"],
+            shape=(3, 4, 5),
+            max_rows=3,
+        )
+
+    assert len(variables) == 2
+    for a in variables:
+        assert isinstance(a, dask.array.Array)
+        assert a.shape == (3, 4, 5)
+        notnull = np.isfinite(a)
+        assert notnull[0, 0, 0]
+        assert notnull[1, 0, 0]
+        assert notnull[2, 0, 0]

From c8ecda8ba0d5610d8f78f289d27f096cfeffda2b Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 19 May 2022 17:14:23 +0200
Subject: [PATCH 20/35] Update disu connectivity tests

---
 imod/tests/test_mf6/test_mf6_disu.py | 26 ++++++++++----------------
 1 file changed, 10 insertions(+), 16 deletions(-)

diff --git a/imod/tests/test_mf6/test_mf6_disu.py b/imod/tests/test_mf6/test_mf6_disu.py
index ed6075540..4b8e1d005 100644
--- a/imod/tests/test_mf6/test_mf6_disu.py
+++ b/imod/tests/test_mf6/test_mf6_disu.py
@@ -29,11 +29,6 @@ def structured_dis():
 
 
 def connectivity_checks(dis):
-    def check_symmetry(i, j, v):
-        order0 = np.lexsort((j, i))
-        order1 = np.lexsort((i, j))
-        assert np.allclose(v[order0], v[order1])
-
     ncell = dis.dataset["node"].size
     i = np.repeat(np.arange(1, ncell + 1), dis.dataset["iac"].values) - 1
     j = dis.dataset["ja"].values - 1
@@ -42,25 +37,20 @@ def check_symmetry(i, j, v):
     hwva = dis.dataset["hwva"].values
     cl12 = dis.dataset["cl12"].values
 
-    node = i == j
-    connection = i != j
+    # The node number itself is marked by a negative number.
+    connection = j > 0
     vertical = (ihc == 0) & connection
 
     assert dis.dataset["iac"].values.sum() == j.size
     assert i.min() == 0
-    assert j.min() == 0
     assert i.max() == ncell - 1
     assert j.max() == ncell - 1
-    assert (diff[node] == 0).all()
     assert (diff[connection] > 0).all()
-    assert (cl12[node] == 0).all()
     assert (cl12[connection] != 0).all()
-    assert (hwva[node] == 0).all()
     assert (hwva[connection] != 0).all()
     assert (diff[vertical] > 1).all()
     assert np.allclose(cl12[vertical], 7.5)
     assert np.allclose(hwva[vertical], 200.0)
-    check_symmetry(i, j, hwva)
 
 
 def test_cell_number():
@@ -97,7 +87,7 @@ def test_structured_connectivity_full():
 
 def test_from_structured(structured_dis):
     idomain, top, bottom = structured_dis
-    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_dis(top, bottom, idomain)
     assert np.allclose(dis.dataset["xorigin"], 40.0)
     assert np.allclose(dis.dataset["yorigin"], 0.0)
     connectivity_checks(dis)
@@ -105,13 +95,13 @@ def test_from_structured(structured_dis):
     # Now disable some cells, create one pass-through
     idomain.values[1, 0, 1] = -1
     idomain.values[:, 0, 0] = 0
-    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_dis(top, bottom, idomain)
     connectivity_checks(dis)
 
 
 def test_render(structured_dis):
     idomain, top, bottom = structured_dis
-    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_dis(top, bottom, idomain)
 
     directory = pathlib.Path("mymodel")
     actual = dis.render(directory, "disu", None, True)
@@ -135,6 +125,8 @@ def test_render(structured_dis):
             open/close mymodel/disu/bot.bin (binary)
           area
             open/close mymodel/disu/area.bin (binary)
+          idomain
+            open/close mymodel/disu/idomain.bin (binary)
         end griddata
 
         begin connectiondata
@@ -150,12 +142,14 @@ def test_render(structured_dis):
             open/close mymodel/disu/hwva.bin (binary)
         end connectiondata"""
     )
+    print(actual)
+    print(expected)
     assert actual == expected
 
 
 def test_write(structured_dis, tmp_path):
     idomain, top, bottom = structured_dis
-    dis = disu.LowLevelUnstructuredDiscretization.from_structured(top, bottom, idomain)
+    dis = disu.LowLevelUnstructuredDiscretization.from_dis(top, bottom, idomain)
     dis.write(tmp_path, "disu", None, True)
 
     assert (tmp_path / "disu.disu").exists

From 836b63d86c6f5771a5ad831e2244a4c216d7c0f7 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 19 May 2022 17:16:58 +0200
Subject: [PATCH 21/35] Use textwrap for test text content

---
 .../test_mf6_input_reading_common.py          | 48 ++++++++--------
 .../test_mf6_input_reading_griddata.py        | 57 ++++++++++---------
 .../test_mf6_input_reading_listinput.py       | 16 +++---
 3 files changed, 63 insertions(+), 58 deletions(-)

diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
index b2b9328b9..2e40d1641 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_common.py
@@ -1,3 +1,5 @@
+import textwrap
+
 import numpy as np
 import pytest
 
@@ -126,17 +128,16 @@ def test_read_external_textfile(tmp_path):
 
 def test_advance_to_header(tmp_path):
     path = tmp_path / "header.txt"
-    content = "\n".join(
-        [
-            "",
-            "begin options",
-            "1",
-            "end options\n",
-            "\n",
-            "begin griddata",
-            "2",
-            "end griddata",
-        ]
+    content = textwrap.dedent(
+        """\
+        begin options
+        1
+        end options
+
+        begin griddata
+        2
+        end griddata
+        """
     )
     with open(path, "w") as f:
         f.write(content)
@@ -213,12 +214,12 @@ def parse(line, _):  # second arg is normally file name
     ]
 
     path2 = tmp_path / "iterable-values-unterminated.txt"
-    content = "\n".join(
-        [
-            "1.0 1 1.0",
-            "1.0 1 1.0",
-            "1.0 1 1.0",
-        ]
+    content = textwrap.dedent(
+        """\
+        1.0 1 1.0
+        1.0 1 1.0
+        1.0 1 1.0
+        """
     )
     with open(path2, "w") as f:
         f.write(content)
@@ -233,15 +234,16 @@ def test_parse_dimension():
 
 def test_advance_to_period(tmp_path):
     path = tmp_path / "period.txt"
-    content = "\n".join(
-        [
-            "" "  begin period 1",
-            "2",
-        ]
+    content = textwrap.dedent(
+        """\
+
+        begin period 1
+        2
+        """
     )
     with open(path, "w") as f:
         f.write(content)
 
     with open(path) as f:
         cm.advance_to_period(f)
-        assert f.readline() == "2"
+        assert f.readline() == "2\n"
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
index f16e668b2..c41eabec7 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_griddata.py
@@ -1,3 +1,5 @@
+import textwrap
+
 import dask.array
 import numpy as np
 import pytest
@@ -7,14 +9,12 @@
 
 def test_advance_to_griddata_section(tmp_path):
     path = tmp_path / "griddata.txt"
-    content = "\n".join(
-        [
-            "",
-            " idomain",
-            "   open/close GWF_1/dis/idomain.dat",
-            " botm layered",
-            "",
-        ]
+    content = textwrap.dedent(
+        """\
+          idomain
+            open/close GWF_1/dis/idomain.dat",
+          botm layered
+        """
     )
 
     with open(path, "w") as f:
@@ -81,16 +81,16 @@ def test_read_array(tmp_path):
     external_path = tmp_path / "external.dat"
     external_binary_path = tmp_path / "external.bin"
 
-    content = "\n".join(
-        [
-            "top",
-            "  constant 200.0",
-            "idomain",
-            "  open/close external.dat",
-            "botm",
-            "  open/close external.bin (binary)",
-            "abc",
-        ]
+    content = textwrap.dedent(
+        """\
+        top
+          constant 200.0
+        idomain
+          open/close external.dat
+        botm
+          open/close external.bin (binary)
+        abc
+        """
     )
     with open(blockfile_path, "w") as f:
         f.write(content)
@@ -126,17 +126,18 @@ def test_read_griddata(tmp_path):
     external_path = tmp_path / "external.dat"
     external_binary_path = tmp_path / "external.bin"
 
-    content = "\n".join(
-        [
-            "top",
-            "  constant 200.0",
-            "idomain",
-            "  open/close external.dat",
-            "botm",
-            "  open/close external.bin (binary)",
-            "end",
-        ]
+    content = textwrap.dedent(
+        """\
+        top
+          constant 200.0
+        idomain
+          open/close external.dat
+        botm
+          open/close external.bin (binary)
+        end
+        """
     )
+
     with open(blockfile_path, "w") as f:
         f.write(content)
 
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
index 473430dcd..b85470b9d 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
@@ -1,3 +1,5 @@
+import textwrap
+
 import dask.array
 import numpy as np
 
@@ -109,13 +111,13 @@ def test_recarr_to_dense__disu():
 def test_read_text_listinput(tmp_path):
     dtype = DIS_DTYPE
     path = tmp_path / "listinput.dat"
-    content = "\n".join(
-        [
-            "# layer row column head conductance",
-            "1 1 1 1.0 10.0",
-            "2 1 1 2.0 10.0",
-            "3 1 1 3.0 10.0",
-        ]
+    content = textwrap.dedent(
+        """\
+        # layer row column head conductance
+        1 1 1 1.0 10.0
+        2 1 1 2.0 10.0
+        3 1 1 3.0 10.0
+        """
     )
 
     with open(path, "w") as f:

From cad182e83a80d7ccccfb14a115475aa62afcf9b2 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 19 May 2022 18:48:39 +0200
Subject: [PATCH 22/35] Forgot to use a tmp_path fixture

---
 .../test_input_reading/test_mf6_input_reading_listinput.py      | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
index b85470b9d..6e0bbf1b3 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
@@ -177,7 +177,7 @@ def test_read_binary_listinput(tmp_path):
 
 def test_read_listinput(tmp_path):
     dtype = DIS_DTYPE
-    path = "package-binary.txt"
+    path = tmp_path / "package-binary.txt"
     content = "\n".join(
         [
             "open/close listinput.bin (binary)",

From 032c88a6faf1665074f20af2e7bc3b77d224e7c2 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 20 May 2022 14:53:46 +0200
Subject: [PATCH 23/35] More unit tests and integration test for mf6 input
 reading

---
 imod/mf6/simulation.py                        |  2 +-
 imod/mf6/sto.py                               |  2 +-
 imod/tests/test_mf6/test_ex01_twri.py         | 25 +++++++++
 .../test_mf6_input_reading.py                 | 52 +++++++++++++++++++
 4 files changed, 79 insertions(+), 2 deletions(-)
 create mode 100644 imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py

diff --git a/imod/mf6/simulation.py b/imod/mf6/simulation.py
index b0be98591..5b86630d2 100644
--- a/imod/mf6/simulation.py
+++ b/imod/mf6/simulation.py
@@ -141,7 +141,7 @@ def open(cls, path: Union[str, pathlib.Path]):
         content = read_sim_namefile(path)
 
         # Get the times from the time discretization file.
-        tdis_path = content["tdis6"]
+        tdis_path = simroot / content["tdis6"]
         tdis = read_tdis(tdis_path)
         globaltimes = tdis_time(tdis)
 
diff --git a/imod/mf6/sto.py b/imod/mf6/sto.py
index f2a898561..ec957afc8 100644
--- a/imod/mf6/sto.py
+++ b/imod/mf6/sto.py
@@ -55,7 +55,7 @@ def open(
             "sy": (np.float64, shape_to_max_rows),
         }
         content = read_sto_blockfile(
-            path=path,
+            path=simroot / path,
             simroot=simroot,
             sections=sections,
             shape=shape,
diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py
index 29cf41c05..1eb3b8434 100644
--- a/imod/tests/test_mf6/test_ex01_twri.py
+++ b/imod/tests/test_mf6/test_ex01_twri.py
@@ -2,6 +2,7 @@
 import sys
 import textwrap
 
+import dask.array
 import numpy as np
 import pytest
 import xarray as xr
@@ -460,3 +461,27 @@ def test_simulation_write_errors(twri_model, tmp_path):
     expected_message = "No sto package found in model GWF_1"
     with pytest.raises(ValueError, match=re.escape(expected_message)):
         simulation.write(modeldir, binary=True)
+
+
+@pytest.mark.usefixtures("transient_twri_model")
+@pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
+def test_simulation_write_and_open(transient_twri_model, tmp_path):
+    simulation = transient_twri_model
+    modeldir = tmp_path / "ex01-twri-transient-binary"
+    simulation.write(modeldir, binary=True)
+
+    back = imod.mf6.Modflow6Simulation.open(modeldir / "mfsim.nam")
+    assert isinstance(back, imod.mf6.Modflow6Simulation)
+
+    gwf = back["gwf_1"]
+    for name in ["chd", "drn", "ic", "npf", "rch", "sto"]:
+        assert name in gwf
+
+    chd = gwf["chd"]
+    assert isinstance(chd, imod.mf6.ConstantHead)
+    assert tuple(chd.dataset["head"].dims) == ("time", "layer", "y", "x")
+    assert isinstance(chd.dataset["head"].data, dask.array.Array)
+
+    head = chd["head"].dropna("layer", how="all").isel(time=0, drop=True).compute()
+    original = simulation["GWF_1"]["chd"]["head"]
+    assert head.equals(original)
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py
new file mode 100644
index 000000000..d57afcd55
--- /dev/null
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py
@@ -0,0 +1,52 @@
+import textwrap
+
+from imod.mf6 import read_input as ri
+
+
+def test_parse_model():
+    ri.parse_model("gwf6 GWF_1.nam GWF_1", "sim.nam") == ["gwf6", "GWF_1_nam", "GWF_1"]
+    with pytest.raises(ValueError, match="ftype, fname and pname expected"):
+        ri.parse_model("gwf6 GWF_1.nam", "sim.nam")
+
+
+def test_parse_exchange():
+    ri.parse_exchange(
+        "GWF6-GWF6 simulation.exg GWF_Model_1 GWF_Model_2", "sim.nam"
+    ) == ["GWF6-GWF6, simulation.exg, GWF_Model_1, GWF_Model_2"]
+    with pytest.raises(ValueError, match="exgtype, exgfile, exgmnamea, exgmnameb"):
+        ri.parse_exchange("GWF6-GWF6 simulation.exg GWF_Model_1", "sim.nam")
+
+
+def test_parse_solutiongroup():
+    ri.parse_solutiongroup("ims6 solver.ims GWF_1", "sim.nam") == [
+        "ims6",
+        "solver.ims",
+        "GWF_1",
+    ]
+    with pytest.raises(ValueError, match="Expected at least three entries"):
+        ri.parse_solutiongroup("ims6 solver.ims", "sim.name")
+
+
+def test_parse_package():
+    ri.parse_package("dis6 dis.dis", "gwf.nam") == ["dis6", "dis.dis", "dis"]
+    ri.parse_package("dis6 dis.dis discretization", "gwf.nam") == [
+        "dis6",
+        "dis.dis",
+        "discretization",
+    ]
+    with pytest.raises(ValueError, match="Expected ftype, fname"):
+        ri.parse_package("dis6", "gwf.nam")
+    with pytest.raises(ValueError, match="Expected ftype, fname"):
+        ri.parse_package("dis6 dis.dis abc def", "gwf.nam")
+
+
+def test_parse_tdis_perioddata():
+    perlen, nstp, tsmult = ri.parse_tdis_perioddata("1.0 1 1.0", "timedis.tdis")
+    assert isinstance(perlen, float)
+    assert perlen == 1.0
+    assert isinstance(nstp, int)
+    assert nstp == 1
+    assert isinstance(tsmult, float)
+    assert tsmult == 1.0
+    with pytest.raises(ValueError, match="perlen, nstp, tsmult expected"):
+        ri.parse_tdis_perioddata("1.0 1.0", "timedis.tdis")

From 76a6b57845265a9a95b38ec5902cb4d2d4eb7087 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 23 May 2022 19:06:59 +0200
Subject: [PATCH 24/35] Add logic for reading model input of classes that are
 DIS, DISV, or DISU specific

---
 imod/mf6/model.py | 52 +++++++++++++++++++++++++++++------------------
 1 file changed, 32 insertions(+), 20 deletions(-)

diff --git a/imod/mf6/model.py b/imod/mf6/model.py
index 096ec2671..2c814bd03 100644
--- a/imod/mf6/model.py
+++ b/imod/mf6/model.py
@@ -30,27 +30,38 @@ class GroundwaterFlowModel(Model):
     _pkg_id = "model"
 
     @staticmethod
-    def _PACKAGE_CLASSES() -> Dict[str, Type]:
+    def _PACKAGE_CLASSES(distype: str) -> Dict[str, Type]:
+        # DIS dependent packages do not have a completely topology description
+        # to a Python equivalent. These are packages such as:
+        #
+        # * Well
+        # * Multi-Aquifer Well
+        # * Horizontal Flow Barrier
+        # * Stream Flow Routing
+        #
+        # Instead, they are returned in a lower level form, directly related to
+        # the MODFLOW6 input.
+        dis_dependent = {
+            "dis": (mf6.WellDisStructured,),
+            "disv": (mf6.WellDisVertices,),
+        }
+
         # mf6.OutputControl is skipped
         # mf6.StorageCoefficient handled by SpecificStorage
-        return {
-            package._pkg_id: package
-            for package in (
-                mf6.ConstantHead,
-                mf6.StructuredDiscretization,
-                mf6.VerticesDiscretization,
-                mf6.Drainage,
-                mf6.Evapotranspiration,
-                mf6.GeneralHeadBoundary,
-                mf6.InitialConditions,
-                mf6.NodePropertyFlow,
-                mf6.Recharge,
-                mf6.River,
-                mf6.SpecificStorage,
-                # mf6.WellDisStructured,
-                # mf6.WellDisVertices,
-            )
-        }
+        packages = dis_dependent[distype] + (
+            mf6.ConstantHead,
+            mf6.StructuredDiscretization,
+            mf6.VerticesDiscretization,
+            mf6.Drainage,
+            mf6.Evapotranspiration,
+            mf6.GeneralHeadBoundary,
+            mf6.InitialConditions,
+            mf6.NodePropertyFlow,
+            mf6.Recharge,
+            mf6.River,
+            mf6.SpecificStorage,
+        )
+        return {package._pkg_id: package for package in packages}
 
     def _initialize_template(self):
         loader = jinja2.PackageLoader("imod", "templates/mf6")
@@ -151,7 +162,6 @@ def open(
 
         # Search for the DIS/DISV/DISU package first. This provides us with
         # the coordinates and dimensions to instantiate the other packages.
-        classes = cls._PACKAGE_CLASSES()
         dis_packages = [
             tup for tup in content["packages"] if tup[0] in ("dis6", "disv6", "disu6")
         ]
@@ -160,6 +170,8 @@ def open(
 
         disftype, disfname, dispname = dis_packages[0]
         diskey = disftype[:-1]
+
+        classes = cls._PACKAGE_CLASSES(diskey)
         package = classes[diskey]
         dis_package = package.open(
             disfname,

From b6ad9d3fb1d8d2144f77e1520c8919d99683e449 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 23 May 2022 19:08:42 +0200
Subject: [PATCH 25/35] Add logic for reading WellDisStructured

* Do not convert to dense (grid) form
* Add optional keyword `sparse_to_dense` in list input reading functions
* Update test to read Well data
---
 imod/mf6/read_input/__init__.py               | 17 ++++++-
 imod/mf6/read_input/list_input.py             | 44 +++++++++++++++----
 imod/mf6/wel.py                               | 26 +++++++++++
 imod/tests/test_mf6/test_ex01_twri.py         |  7 +--
 .../test_mf6_input_reading.py                 |  2 +-
 .../test_mf6_input_reading_listinput.py       | 23 ++++++++++
 6 files changed, 104 insertions(+), 15 deletions(-)

diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index e87b75d5b..f788c2e45 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -274,6 +274,7 @@ def read_package_periods(
     fields: Tuple[str],
     max_rows: int,
     shape: Tuple[int],
+    sparse_to_dense: bool = True,
 ) -> Tuple[List[int], Dict[str, dask.array.Array]]:
     """
     Read blockfile periods section of a "standard" MODFLOW6 boundary
@@ -308,6 +309,11 @@ def read_package_periods(
     dask_lists = defaultdict(list)
     key = advance_to_period(f)
 
+    if sparse_to_dense:
+        variables = fields
+    else:
+        variables = index_columns + fields
+
     while key is not None:
         # Read the recarrays, already separated into dense arrays.
         variable_values = read_listinput(
@@ -318,11 +324,12 @@ def read_package_periods(
             index_columns=index_columns,
             max_rows=max_rows,
             shape=shape,
+            sparse_to_dense=sparse_to_dense,
         )
 
         # Group them by field (e.g. cond, head, etc.)
-        for field, values in zip(fields, variable_values):
-            dask_lists[field].append(values)
+        for var, values in zip(variables, variable_values):
+            dask_lists[var].append(values)
 
         # Store number and advance to next period
         periods.append(key)
@@ -341,6 +348,7 @@ def read_boundary_blockfile(
     simroot: Path,
     fields: Tuple[str],
     shape: Tuple[int],
+    sparse_to_dense: bool = True,
 ) -> Dict[str, Any]:
     """
     Read blockfile of a "standard" MODFLOW6 boundary condition package: RIV,
@@ -358,6 +366,10 @@ def read_boundary_blockfile(
         The fields (generally np.float64) of the dtype.
     shape: Tuple[int]
         DIS: 3-sized, DISV: 2-sized, DISU: 1-sized.
+    sparse_to_dense: bool, default: True
+        Whether to convert "sparse" COO (cell id) data into "dense" grid form
+        (with implicit location). False for packages such as Well, which are
+        not usually in grid form.
 
     Returns
     -------
@@ -408,6 +420,7 @@ def read_boundary_blockfile(
             fields=fields,
             max_rows=dimensions["maxbound"],
             shape=shape,
+            sparse_to_dense=sparse_to_dense,
         )
 
     return content
diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
index 1307e6eb1..7222ddbeb 100644
--- a/imod/mf6/read_input/list_input.py
+++ b/imod/mf6/read_input/list_input.py
@@ -11,6 +11,16 @@
 from .common import read_external_binaryfile, split_line, strip_line
 
 
+def field_values(
+    recarr: np.ndarray,
+    fields: Tuple[str],
+):
+    """
+    Return the record array columns as a list of separate arrays.
+    """
+    return [recarr[field] for field in fields]
+
+
 def recarr_to_dense(
     recarr: np.ndarray,
     index_columns: Tuple[str],
@@ -24,9 +34,9 @@ def recarr_to_dense(
     # MODFLOW6 is 1-based, Python is 0-based
     indices = [recarr[column] - 1 for column in index_columns]
     variables = []
-    for field in fields:
+    for column in field_values(recarr, fields):
         data = np.full(shape, np.nan)
-        data[indices] = recarr[field]
+        data[indices] = column
         variables.append(data)
     return variables
 
@@ -59,10 +69,14 @@ def read_internal_listinput(
     fields: Tuple[str],
     shape: Tuple[int],
     max_rows: int,
+    sparse_to_dense: bool,
 ) -> List[np.ndarray]:
     # recarr = read_internal(f, max_rows, dtype)
     recarr = read_text_listinput(f, dtype, max_rows)
-    return recarr_to_dense(recarr, index_columns, fields, shape)
+    if sparse_to_dense:
+        return recarr_to_dense(recarr, index_columns, fields, shape)
+    else:
+        return field_values(recarr, index_columns + fields)
 
 
 def read_external_listinput(
@@ -73,6 +87,7 @@ def read_external_listinput(
     shape: Tuple[int],
     binary: bool,
     max_rows: int,
+    sparse_to_dense: bool,
 ) -> List[np.ndarray]:
     """
     Read external list input, separate and reshape to a dense array form.
@@ -81,7 +96,10 @@ def read_external_listinput(
         recarr = read_external_binaryfile(path, dtype, max_rows)
     else:
         recarr = read_text_listinput(path, dtype, max_rows)
-    return recarr_to_dense(recarr, index_columns, fields, shape)
+    if sparse_to_dense:
+        return recarr_to_dense(recarr, index_columns, fields, shape)
+    else:
+        return field_values(recarr, index_columns + fields)
 
 
 def read_listinput(
@@ -92,6 +110,7 @@ def read_listinput(
     fields: Tuple[str],
     shape: Tuple[int],
     max_rows: int,
+    sparse_to_dense: bool = True,
 ) -> List[dask.array.Array]:
     """
     MODFLOW6 list input reading functionality.
@@ -111,6 +130,7 @@ def read_listinput(
     fields: Tuple[str]
     shape: Tuple[int]
     max_rows: int
+    sparse_to_dense: bool
 
     Returns
     -------
@@ -127,6 +147,13 @@ def read_listinput(
     separated = split_line(stripped)
     first = separated[0]
 
+    nout = len(fields)
+    fieldtypes = [dtype.fields[field][0] for field in fields]
+    if not sparse_to_dense:
+        shape = (max_rows,)
+        nout += len(index_columns)
+        fieldtypes = [dtype.fields[field][0] for field in index_columns] + fieldtypes
+
     if first == "open/close":
         fname = separated[1]
         path = simroot / fname
@@ -135,12 +162,9 @@ def read_listinput(
         if binary and "boundnames" in dtype.fields:
             raise ValueError("(BINARY) does not support BOUNDNAMES")
 
-        nout = len(fields)
         x = dask.delayed(read_external_listinput, nout=nout)(
-            path, dtype, index_columns, fields, shape, binary, max_rows
+            path, dtype, index_columns, fields, shape, binary, max_rows, sparse_to_dense
         )
-
-        fieldtypes = [dtype.fields[field][0] for field in fields]
         variable_values = [
             dask.array.from_delayed(a, shape=shape, dtype=dt)
             for a, dt in zip(x, fieldtypes)
@@ -148,7 +172,9 @@ def read_listinput(
     else:
         # Move file position back one line, and try reading internal values.
         f.seek(position)
-        x = read_internal_listinput(f, dtype, index_columns, fields, shape, max_rows)
+        x = read_internal_listinput(
+            f, dtype, index_columns, fields, shape, max_rows, sparse_to_dense
+        )
         variable_values = [dask.array.from_array(a) for a in x]
 
     return variable_values
diff --git a/imod/mf6/wel.py b/imod/mf6/wel.py
index b0c77af3d..551436a6f 100644
--- a/imod/mf6/wel.py
+++ b/imod/mf6/wel.py
@@ -5,6 +5,8 @@
 
 from imod.mf6.pkgbase import BoundaryCondition, VariableMetaData
 
+from .read_input import read_boundary_blockfile
+
 
 def assign_dims(arg) -> Dict:
     is_da = isinstance(arg, xr.DataArray)
@@ -109,6 +111,30 @@ def to_sparse(self, arrdict, layer):
             recarr[key] = arr
         return recarr
 
+    @classmethod
+    def open(cls, path, simroot, shape, coords, dims, globaltimes):
+        # read_boundary_blockfile calls:
+        # read_package_periods, which calls:
+        # read_listinput, which transforms recarrays to dense arrays by default.
+        content = read_boundary_blockfile(
+            simroot / path,
+            simroot,
+            fields=("rate",),
+            shape=shape,
+            sparse_to_dense=False,
+        )
+
+        period_index = content.pop("period_index")
+        period_data = content.pop("period_data")
+        coords = {"time": globaltimes[period_index]}
+        dims = ("time", "index")
+
+        for field, data in period_data.items():
+            content[field] = xr.DataArray(data, coords, dims)
+
+        filtered_content = cls.filter_and_rename(content)
+        return cls(**filtered_content)
+
 
 class WellDisVertices(BoundaryCondition):
     """
diff --git a/imod/tests/test_mf6/test_ex01_twri.py b/imod/tests/test_mf6/test_ex01_twri.py
index 1eb3b8434..01eac3637 100644
--- a/imod/tests/test_mf6/test_ex01_twri.py
+++ b/imod/tests/test_mf6/test_ex01_twri.py
@@ -463,18 +463,19 @@ def test_simulation_write_errors(twri_model, tmp_path):
         simulation.write(modeldir, binary=True)
 
 
+@pytest.mark.parametrize("binary", [True, False])
 @pytest.mark.usefixtures("transient_twri_model")
 @pytest.mark.skipif(sys.version_info < (3, 7), reason="capture_output added in 3.7")
-def test_simulation_write_and_open(transient_twri_model, tmp_path):
+def test_simulation_write_and_open(transient_twri_model, tmp_path, binary):
     simulation = transient_twri_model
     modeldir = tmp_path / "ex01-twri-transient-binary"
-    simulation.write(modeldir, binary=True)
+    simulation.write(modeldir, binary=binary)
 
     back = imod.mf6.Modflow6Simulation.open(modeldir / "mfsim.nam")
     assert isinstance(back, imod.mf6.Modflow6Simulation)
 
     gwf = back["gwf_1"]
-    for name in ["chd", "drn", "ic", "npf", "rch", "sto"]:
+    for name in ["chd", "drn", "ic", "npf", "rch", "sto", "wel"]:
         assert name in gwf
 
     chd = gwf["chd"]
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py
index d57afcd55..88df4f67e 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading.py
@@ -1,4 +1,4 @@
-import textwrap
+import pytest
 
 from imod.mf6 import read_input as ri
 
diff --git a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
index 6e0bbf1b3..87925119a 100644
--- a/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
+++ b/imod/tests/test_mf6/test_input_reading/test_mf6_input_reading_listinput.py
@@ -132,6 +132,7 @@ def test_read_text_listinput(tmp_path):
             fields=["head", "conductance"],
             max_rows=3,
             shape=(3, 4, 5),
+            sparse_to_dense=True,
         )
     assert isinstance(variables, list)
     assert len(variables) == 2
@@ -144,6 +145,7 @@ def test_read_text_listinput(tmp_path):
         shape=(3, 4, 5),
         binary=False,
         max_rows=3,
+        sparse_to_dense=True,
     )
     assert isinstance(variables, list)
     assert len(variables) == 2
@@ -170,6 +172,7 @@ def test_read_binary_listinput(tmp_path):
         shape=(3, 4, 5),
         binary=True,
         max_rows=3,
+        sparse_to_dense=True,
     )
     assert isinstance(variables, list)
     assert len(variables) == 2
@@ -217,3 +220,23 @@ def test_read_listinput(tmp_path):
         assert notnull[0, 0, 0]
         assert notnull[1, 0, 0]
         assert notnull[2, 0, 0]
+
+    # Test another round, this time without converting from COO to a dense
+    # form.
+    with open(path) as f:
+        variables = li.read_listinput(
+            f,
+            tmp_path,
+            dtype,
+            index_columns=["layer", "row", "column"],
+            fields=["head", "conductance"],
+            shape=(3, 4, 5),
+            max_rows=3,
+            sparse_to_dense=False,
+        )
+
+    assert len(variables) == 5
+    for a in variables:
+        assert isinstance(a, dask.array.Array)
+        assert a.shape == (3,)
+        assert np.isfinite(a).all()

From 355531a0b7aa54a4ae659ae4862534d815bbfa44 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 25 May 2022 15:38:51 +0200
Subject: [PATCH 26/35] Move dupuit example to a separate repo

---
 examples/mf6/dupuit_forchheimer_hack.py | 439 ------------------------
 1 file changed, 439 deletions(-)
 delete mode 100644 examples/mf6/dupuit_forchheimer_hack.py

diff --git a/examples/mf6/dupuit_forchheimer_hack.py b/examples/mf6/dupuit_forchheimer_hack.py
deleted file mode 100644
index d12e0648e..000000000
--- a/examples/mf6/dupuit_forchheimer_hack.py
+++ /dev/null
@@ -1,439 +0,0 @@
-"""
-Dupuit-Forchheimer "hack"
-=========================
-
-In contrast to earlier versions of MODFLOW, MODFLOW6 only officially supports
-fully 3D simulations. This means that the number of layers doubles minus one:
-before, heads were only computed for aquifers while in MODFLOW6 heads are
-computed for both aquifers and aquitards.
-
-This has some downsides: the required amount of RAM increases, and computation
-times mildly increase (by roughly 70%). However, MODFLOW6's logic can be
-modified relatively easy through its API, provided by ``xmipy`` in Python.
-This requires a script utilizing ``xmipy``, and customized MODFLOW6 input.
-This example demonstrates both.
-
-XMI script
-==========
-
-In a nutshell, this hack consists of:
-
-* Write MODFLOW6 input with reduced layers, representing the aquifers
-  exclusively.
-* As we cannot rely on MODFLOW6 to use the aquitard thickness to compute
-  aquitard resistance, we assign it ourselves to the K33 entry of the
-  NodePropertyFlow package.
-* Initialize MODFLOW6, fetch the K33 values from memory, and assign these
-  to the saturated conductance (CONDSAT).
-* Run the MODFLOW6 simulation to completion.
-
-We'll start by implementing a DupuitForchheimerSimulation class.
-"""
-#%%
-
-import datetime
-
-import numpy as np
-import pandas as pd
-import xarray as xr
-from xmipy import XmiWrapper
-
-import imod
-
-# %%
-# Some type hinting.
-
-FloatArray = np.ndarray
-IntArray = np.ndarray
-BoolArray = np.ndarray
-
-# %%
-# Simulation class
-# ----------------
-#
-# We'll start by defining an ordinary Simulation class. This class behaves the
-# same as a regular MODFLOW6 simulation (i.e. when running ``mf6.exe```).
-
-
-class Simulation:
-    """
-    Run all stress periods in a simulation.
-    """
-
-    def __init__(self, wdir: str, name: str):
-        self.modelname = name
-        self.mf6 = XmiWrapper(lib_path="libmf6.dll", working_directory=wdir)
-        self.mf6.initialize()
-        self.max_iter = self.mf6.get_value_ptr("SLN_1/MXITER")[0]
-        shape = np.zeros(1, dtype=np.int32)
-        self.ncell = self.mf6.get_grid_shape(1, shape)[0]
-
-    def __enter__(self):
-        return self
-
-    def __exit__(self, type, value, traceback):
-        self.finalize()
-
-    def do_iter(self, sol_id: int) -> bool:
-        """Execute a single iteration"""
-        has_converged = self.mf6.solve(sol_id)
-        return has_converged
-
-    def update(self):
-        # We cannot set the timestep (yet) in Modflow
-        # -> set to the (dummy) value 0.0 for now
-        self.mf6.prepare_time_step(0.0)
-        self.mf6.prepare_solve(1)
-        # Convergence loop
-        self.mf6.prepare_solve(1)
-        for kiter in range(1, self.max_iter + 1):
-            has_converged = self.do_iter(1)
-            if has_converged:
-                print(f"MF6 converged in {kiter} iterations")
-                break
-        self.mf6.finalize_solve(1)
-        # Finish timestep
-        self.mf6.finalize_time_step()
-        current_time = self.mf6.get_current_time()
-        return current_time
-
-    def get_times(self):
-        """Return times"""
-        return (
-            self.mf6.get_start_time(),
-            self.mf6.get_current_time(),
-            self.mf6.get_end_time(),
-        )
-
-    def run(self):
-        start = datetime.datetime.now()
-
-        _, current_time, end_time = self.get_times()
-        while current_time < end_time:
-            current_time = self.update()
-
-        stop = datetime.datetime.now()
-        print(
-            f"Elapsed run time: {stop - start} (hours: minutes: seconds)."
-            f"Simulation terminated normally."
-        )
-
-    def finalize(self):
-        self.mf6.finalize()
-
-
-# %%
-# Dupuit Forchheimer Simulation
-# -----------------------------
-#
-# Next, we'll inherit from the basic simulation and add some additional methods
-# to extract the K33 values, and use it to set the CONDSAT array of the model.
-
-
-class DupuitForchheimerSimulation(Simulation):
-    def __init__(self, wdir: str, name: str):
-        super().__init__(wdir, name)
-        self.set_resistance()
-
-    def conductance_index(
-        self,
-        vertical: BoolArray,
-    ):
-        """
-        We've looked at the upper diagonal half of the coefficient matrix.
-
-        While this means that j is always larger than i, it does not mean that cell
-        i is always overlying cell j; the cell numbering may be arbitrary.
-
-        In our convention, the k33 values is interpreted as the resistance to the
-        cell BELOW it. Should now a case arise when j > i, but with cell j overlying
-        cell i, we should use the k33 value of cell j to set the conductance.
-
-        # TODO: reduced numbers?
-        """
-        mf6 = self.mf6
-        modelname = self.modelname
-        top = mf6.get_value_ptr(f"{modelname}/DIS/TOP")
-        bottom = mf6.get_value_ptr(f"{modelname}/DIS/BOT")
-        # Collect cell-to-cell connections
-        # Python is 0-based, Fortran is 1-based
-        # TODO: grab IAUSR and JAUSR if NODEREDUCED?
-        ia = mf6.get_value_ptr(f"{modelname}/CON/IA") - 1
-        ja = mf6.get_value_ptr(f"{modelname}/CON/JA") - 1
-        # Convert compressed sparse row (CSR) to row(i) and column(j) numbers:
-        n = np.diff(ia)
-        i = np.repeat(np.arange(self.ncell), n)
-        j = ja
-        # Get the upper diagonal, and only the vertical entries.
-        upper = j > i
-        i = i[upper][vertical]
-        j = j[upper][vertical]
-        # Now find out which cell is on top, i or j.
-        top_i = top[i]
-        bottom_j = bottom[j]
-        take_i = top_i >= bottom_j
-        take_j = ~take_i
-        # Create the index with the appropriate values of i and j.
-        index = np.empty(i.size, dtype=np.int32)
-        index[take_i] = i[take_i]
-        index[take_j] = j[take_j]
-        return index
-
-    def set_resistance(self):
-        mf6 = self.mf6
-        modelname = self.modelname
-        # Grab views on the MODFLOW6 memory:
-        area = mf6.get_value_ptr(f"{modelname}/DIS/AREA")
-        k33 = mf6.get_value_ptr(f"{modelname}/NPF/K33")
-        condsat = mf6.get_value_ptr(f"{modelname}/NPF/CONDSAT")
-        ihc = mf6.get_value_ptr(f"{modelname}/CON/IHC")
-        vertical = ihc == 0
-        new_cond = area / k33
-        cell_i = self.conductance_index(vertical)
-        condsat[vertical] = new_cond[cell_i]
-
-
-# %%
-# Creating input
-# ==============
-#
-# As mentioned, the primary thing is to write the resistance of the aquitard
-# layers into K33. Unfortunately, there is one more hurdle: in MODFLOW6, all
-# layers are assumed to be contiguous in depth for DIS and DISV discretization
-# packages. This means that for a layer, the bottom of the overlying layer is
-# equal its top. DISU is an exception, and allows separate top and bottom
-# values. Consequently, we will write all the model as if it is a DISU model.
-#
-# We can do with relative easy, by creating a
-# ``LowLevelUnstructuredDiscretization`` instance via its ``from_dis`` method.
-# Additionally, all packages features a ``to_disu`` method. Otherwise, we
-# assemble and write the model as usual.
-#
-# However, let's start with a small fully 3D benchmark. The model features:
-#
-# * recharge across the entire top layer;
-# * two aquifers of 50 m separated by an aquitard of 10 m;
-# * a drain in the center.
-
-nlay = 3
-nrow = 150
-ncol = 151
-shape = (nlay, nrow, ncol)
-
-dx = 10.0
-dy = -10.0
-xmin = 0.0
-xmax = dx * ncol
-ymin = 0.0
-ymax = abs(dy) * nrow
-dims = ("layer", "y", "x")
-
-layer = np.array([1, 2, 3])
-y = np.arange(ymax, ymin, dy) + 0.5 * dy
-x = np.arange(xmin, xmax, dx) + 0.5 * dx
-coords = {"layer": layer, "y": y, "x": x}
-
-idomain = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
-bottom = xr.DataArray([40.0, 30.0, 0.0], {"layer": layer}, ("layer",))
-top = xr.DataArray([50.0])
-icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",))
-recharge = xr.full_like(idomain.sel(layer=1), 0.001)
-
-# We'll assume a resistance of 10 days, for a ditch of 2 m wide.
-conductance = xr.full_like(idomain.sel(layer=1), np.nan)
-elevation = xr.full_like(idomain.sel(layer=1), np.nan)
-conductance[:, 75] = 20.0
-elevation[:, 75] = 42.0
-
-# %%
-# For comparison, we set horizontal conductivity of the aquitard to a tiny value.
-
-k = xr.DataArray([10.0, 1.0e-6, 10.0], {"layer": layer}, ("layer",))
-k33 = xr.DataArray([10.0, 0.1, 10.0], {"layer": layer}, ("layer",))
-
-gwf_model = imod.mf6.GroundwaterFlowModel()
-gwf_model["dis"] = imod.mf6.StructuredDiscretization(
-    top=top, bottom=bottom, idomain=idomain
-)
-gwf_model["drn"] = imod.mf6.Drainage(
-    elevation=elevation,
-    conductance=conductance,
-    print_input=True,
-    print_flows=True,
-    save_flows=True,
-)
-gwf_model["ic"] = imod.mf6.InitialConditions(head=45.0)
-gwf_model["npf"] = imod.mf6.NodePropertyFlow(
-    icelltype=icelltype,
-    k=k,
-    k33=k33,
-)
-gwf_model["sto"] = imod.mf6.SpecificStorage(
-    specific_storage=1.0e-5,
-    specific_yield=0.15,
-    transient=False,
-    convertible=0,
-)
-gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all")
-gwf_model["rch"] = imod.mf6.Recharge(recharge)
-
-# Attach it to a simulation
-simulation = imod.mf6.Modflow6Simulation("ex01-twri")
-simulation["GWF_1"] = gwf_model
-# Define solver settings
-simulation["solver"] = imod.mf6.Solution(
-    print_option="summary",
-    csv_output=False,
-    no_ptc=True,
-    outer_dvclose=1.0e-4,
-    outer_maximum=500,
-    under_relaxation=None,
-    inner_dvclose=1.0e-4,
-    inner_rclose=0.001,
-    inner_maximum=100,
-    linear_acceleration="cg",
-    scaling_method=None,
-    reordering_method=None,
-    relaxation_factor=0.97,
-)
-simulation.time_discretization(["2020-01-01", "2020-01-02"])
-
-# %%
-
-simulation.write("model3d")
-# %%
-# Let's run the simulation with our class defined above.
-
-with Simulation("model3d", "model3d") as xmi_simulation:
-    xmi_simulation.run()
-
-# %%
-
-head3d = imod.mf6.open_hds_like("model3d/GWF_1/GWF_1.hds", idomain)
-
-# %%
-head3d.isel(time=0, layer=0).plot.contourf()
-
-# %%
-
-headdiff = head3d.isel(time=0, layer=0) - head3d.isel(time=0, layer=1)
-headdiff.plot.imshow()
-
-# %%
-# We can now create the reduced model.
-
-nlay = 2
-nrow = 150
-ncol = 151
-
-layer = np.array([1, 2])
-y = np.arange(ymax, ymin, dy) + 0.5 * dy
-x = np.arange(xmin, xmax, dx) + 0.5 * dx
-coords = {"layer": layer, "y": y, "x": x}
-
-shape = (nlay, nrow, ncol)
-idomain = xr.DataArray(np.ones(shape), coords=coords, dims=dims)
-bottom = idomain * xr.DataArray([40.0, 0.0], {"layer": layer}, ("layer",))
-top = idomain * xr.DataArray([50.0, 30.0], {"layer": layer}, ("layer",))
-icelltype = idomain * xr.DataArray([1, 0], {"layer": layer}, ("layer",))
-recharge = xr.full_like(idomain.sel(layer=1), 0.001)
-
-# We'll assume a resistance of 10 days, for a ditch of 2 m wide.
-conductance = xr.full_like(idomain.sel(layer=1), np.nan)
-elevation = xr.full_like(idomain.sel(layer=1), np.nan)
-conductance[:, 75] = 20.0
-elevation[:, 75] = 42.0
-
-# %%
-# For comparison, we set horizontal conductivity of the aquitard to a tiny value.
-
-k = idomain * xr.DataArray([10.0, 10.0], {"layer": layer}, ("layer",))
-k33 = idomain * xr.DataArray([100.0, 100.0], {"layer": layer}, ("layer",))
-
-gwf_model = imod.mf6.GroundwaterFlowModel()
-gwf_model["disu"] = imod.mf6.LowLevelUnstructuredDiscretization.from_dis(
-    top=top, bottom=bottom, idomain=idomain
-)
-# %%
-gwf_model["drn"] = imod.mf6.Drainage(
-    elevation=elevation,
-    conductance=conductance,
-    print_input=True,
-    print_flows=True,
-    save_flows=True,
-).to_disu()
-gwf_model["ic"] = imod.mf6.InitialConditions(head=45.0)
-gwf_model["npf"] = imod.mf6.NodePropertyFlow(
-    icelltype=icelltype,
-    k=k,
-    k33=k33,
-).to_disu()
-gwf_model["sto"] = imod.mf6.SpecificStorage(
-    specific_storage=1.0e-5,
-    specific_yield=0.15,
-    transient=False,
-    convertible=0,
-)
-gwf_model["oc"] = imod.mf6.OutputControl(save_head="all", save_budget="all")
-gwf_model["rch"] = imod.mf6.Recharge(recharge).to_disu()
-
-# Attach it to a simulation
-simulation = imod.mf6.Modflow6Simulation("dupuit")
-simulation["GWF_1"] = gwf_model
-# Define solver settings
-simulation["solver"] = imod.mf6.Solution(
-    print_option="summary",
-    csv_output=False,
-    no_ptc=True,
-    outer_dvclose=1.0e-4,
-    outer_maximum=500,
-    under_relaxation=None,
-    inner_dvclose=1.0e-4,
-    inner_rclose=0.001,
-    inner_maximum=100,
-    linear_acceleration="cg",
-    scaling_method=None,
-    reordering_method=None,
-    relaxation_factor=0.97,
-)
-simulation.time_discretization(["2020-01-01", "2020-01-02"])
-# %%
-
-simulation.write("dupuit")
-
-# %%
-
-with DupuitForchheimerSimulation("dupuit", "GWF_1") as xmi_sim:
-    xmi_sim.run()
-
-# %%
-
-ncell = 45300
-path = "dupuit/GWF_1/GWF_1.hds"
-disu_coords = {"node": np.arange(1, ncell + 1)}
-d = {"ncell": ncell, "name": "head", "coords": disu_coords}
-head = imod.mf6.out.disu.open_hds(path, d, True)
-
-# %%
-from typing import Any, Dict
-
-
-def unstack(da: xr.DataArray, dim: str, coords: Dict[str, Any]):
-    """
-    Unstack existing dimension into multiple new dimensions with new
-    coordinates.
-    """
-    new = da.copy()
-    new.coords[dim] = pd.MultiIndex.from_product(
-        [v.values for v in coords.values()], names=list(coords.keys())
-    )
-    return new.unstack(dim)
-
-
-head_df = unstack(head, "node", dict(idomain.coords))
-
-# %%
-
-delta = head3d.isel(time=0, layer=0) - head_df.isel(time=0, layer=0)
-
-# %%

From b7c7b3c268ab841e07a0dfbbf7dffdb65d752a60 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 25 May 2022 16:34:48 +0200
Subject: [PATCH 27/35] Add to_disu tests

---
 imod/mf6/pkgbase.py                 | 16 +++++++++---
 imod/tests/test_mf6/test_mf6_drn.py | 20 +++++++++++++++
 imod/tests/test_mf6/test_mf6_npf.py | 40 +++++++++++++++++++++++++++++
 3 files changed, 73 insertions(+), 3 deletions(-)

diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index 58559e780..1e5b41a6d 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -327,7 +327,9 @@ def render(self, directory, pkgname, globaltimes, binary):
     @staticmethod
     def _is_xy_data(obj):
         if isinstance(obj, (xr.DataArray, xr.Dataset)):
-            # TODO
+            # "nja" is not really xy_data: arguably the method name should be
+            # changed. This method is exclusively used whether the contents of
+            # the object should be written to an external file.
             xy = ("x" in obj.dims and "y" in obj.dims) or (
                 "node" in obj.dims or "nja" in obj.dims
             )
@@ -490,6 +492,8 @@ def write_netcdf(self, directory, pkgname, aggregate_layers=False):
 
     def _dis_to_disu(self, cell_ids=None):
         structured = self.dataset
+        if cell_ids is not None and not np.issubdtype(cell_ids.dtype, np.integer):
+            raise TypeError(f"cell_ids should be integer, received: {cell_ids.dtype}")
         if "layer" not in structured.coords:
             raise ValueError("layer coordinate is required")
         if "layer" not in structured.dims:
@@ -503,14 +507,20 @@ def _dis_to_disu(self, cell_ids=None):
         index = np.add.outer(offset, np.arange(ncell_per_layer)).ravel()
 
         if cell_ids is not None:
+            # node has the shape of index, i.e. one value for every cell in DIS
+            # form, including no data padding, which are values of -1.
             node = cell_ids[index]
+            # Create a subselection without the inactive cells.
+            active = [node != -1]
+            dataset = dataset.isel(node=index[active])
+            dataset = dataset.assign_coords(node=node[active])
         else:
-            node = index
+            dataset = dataset.assign_coords(node=index)
 
-        dataset = dataset.assign_coords(node=node)
         return self.__class__(**dataset)
 
     def to_disu(self, cell_ids=None):
+
         """
         Parameters
         ----------
diff --git a/imod/tests/test_mf6/test_mf6_drn.py b/imod/tests/test_mf6/test_mf6_drn.py
index 206108f09..f380833e3 100644
--- a/imod/tests/test_mf6/test_mf6_drn.py
+++ b/imod/tests/test_mf6/test_mf6_drn.py
@@ -84,3 +84,23 @@ def test_3d_singelayer():
     arrdict = drn._ds_to_arrdict(bin_ds)
     sparse_data = drn.to_sparse(arrdict, layer)
     assert isinstance(sparse_data, np.ndarray)
+
+
+def test_to_disu(drainage):
+    drn = drainage
+    disu_drn = drn.to_disu()
+    assert tuple(disu_drn["conductance"].dims) == ("node",)
+    assert np.array_equal(disu_drn.dataset["node"], np.arange(75))
+
+    drn["layer"] = [1, 3, 5]
+    disu_drn = drn.to_disu()
+    expected = np.concatenate([np.arange(25), np.arange(50, 75), np.arange(100, 125)])
+    assert np.array_equal(disu_drn.dataset["node"], expected)
+
+    onelayer = imod.mf6.Drainage(**drainage.dataset.isel(layer=0))
+    disu_drn = onelayer.to_disu()
+    assert np.array_equal(disu_drn.dataset["node"], np.arange(25))
+
+    nolayer = imod.mf6.Drainage(**drainage.dataset.isel(layer=0, drop=True))
+    with pytest.raises(ValueError, match="layer coordinate is required"):
+        nolayer.to_disu()
diff --git a/imod/tests/test_mf6/test_mf6_npf.py b/imod/tests/test_mf6/test_mf6_npf.py
index 735f16897..607807017 100644
--- a/imod/tests/test_mf6/test_mf6_npf.py
+++ b/imod/tests/test_mf6/test_mf6_npf.py
@@ -68,3 +68,43 @@ def test_wrong_dtype():
             perched=True,
             save_flows=True,
         )
+
+
+def test_to_disu():
+    layer = [1, 2, 3]
+    template = xr.full_like(
+        imod.util.empty_3d(10.0, 0.0, 100.0, 10.0, 0.0, 100.0, layer), 1, dtype=np.int32
+    )
+    icelltype = xr.DataArray([1, 0, 0], {"layer": layer}, ("layer",)) * template
+    k = xr.DataArray([1.0e-3, 1.0e-4, 2.0e-4], {"layer": layer}, ("layer",)) * template
+    k33 = (
+        xr.DataArray([2.0e-8, 2.0e-8, 2.0e-8], {"layer": layer}, ("layer",)) * template
+    )
+
+    npf = imod.mf6.NodePropertyFlow(
+        icelltype=icelltype,
+        k=k,
+        k33=k33,
+        variable_vertical_conductance=True,
+        dewatered=True,
+        perched=True,
+        save_flows=True,
+    )
+
+    disu_npf = npf.to_disu()
+    assert np.array_equal(disu_npf.dataset["node"], np.arange(300))
+
+    # Test again, with partially active domain
+    active = xr.full_like(template, True, dtype=bool)
+    active[:, :, 0] = False
+    n_active = active.sum()
+    cell_ids = np.full((3, 10, 10), -1)
+    cell_ids[active.values] = np.arange(n_active)
+    cell_ids = cell_ids.ravel()
+
+    disu_npf = npf.to_disu(cell_ids)
+    assert np.array_equal(disu_npf.dataset["node"], np.arange(270))
+
+    with pytest.raises(TypeError, match="cell_ids should be integer"):
+        cell_ids = np.arange(3 * 5 * 5, dtype=np.float64)
+        npf.to_disu(cell_ids)

From 35976e329bfa5123b31ef3dfb60999807f4854b0 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 25 May 2022 16:43:03 +0200
Subject: [PATCH 28/35] Remove unused bas module from mf6

---
 imod/mf6/bas.py | 215 ------------------------------------------------
 1 file changed, 215 deletions(-)
 delete mode 100644 imod/mf6/bas.py

diff --git a/imod/mf6/bas.py b/imod/mf6/bas.py
deleted file mode 100644
index 265a5d44a..000000000
--- a/imod/mf6/bas.py
+++ /dev/null
@@ -1,215 +0,0 @@
-import jinja2
-import numpy as np
-import scipy.ndimage.morphology
-import xarray as xr
-
-from imod import util
-from imod.wq.pkgbase import Package
-
-
-class BasicFlow(Package):
-    """
-    The Basic package is used to specify certain data used in all models.
-    These include:
-
-    1. the locations of acitve, inactive, and specified head in cells,
-    2. the head stored in inactive cells,
-    3. the initial head in all cells, and
-    4. the top and bottom of the aquifer
-
-    The number of layers (NLAY) is automatically calculated using the IBOUND.
-    Thickness is calculated using the specified tops en bottoms.
-    The Basic package input file is required in all models.
-
-    Parameters
-    ----------
-    ibound: xr.DataArray of integers
-        is the boundary variable.
-        If IBOUND(J,I,K) < 0, cell J,I,K has a constant head.
-        If IBOUND(J,I,K) = 0, cell J,I,K is inactive.
-        If IBOUND(J,I,K) > 0, cell J,I,K is active.
-    top: float or xr.DataArray of floats
-        is the top elevation of layer 1. For the common situation in which the
-        top layer represents a water-table aquifer, it may be reasonable to set
-        `top` equal to land-surface elevation.
-    bottom: xr.DataArray of floats
-        is the bottom elevation of model layers or Quasi-3d confining beds. The
-        DataArray should at least include the `layer` dimension.
-    starting_head: float or xr.DataArray of floats
-        is initial (starting) head—that is, head at the beginning of the
-        simulation (STRT). starting_head must be specified for all simulations,
-        including steady-state simulations. One value is read for every model
-        cell. Usually, these values are read a layer at a time.
-    inactive_head: float, optional
-        is the value of head to be assigned to all inactive (no flow) cells
-        (IBOUND = 0) throughout the simulation (HNOFLO). Because head at
-        inactive cells is unused in model calculations, this does not affect
-        model results but serves to identify inactive cells when head is
-        printed. This value is also used as drawdown at inactive cells if the
-        drawdown option is used. Even if the user does not anticipate having
-        inactive cells, a value for inactive_head must be entered. Default
-        value is 1.0e30.
-    """
-
-    _pkg_id = "bas6"
-    _template = jinja2.Template(
-        "[bas6]\n"
-        "    {%- for layer, value in ibound.items() %}\n"
-        "    ibound_l{{layer}} = {{value}}\n"
-        "    {%- endfor %}\n"
-        "    hnoflo = {{inactive_head}}\n"
-        "    {%- for layer, value in starting_head.items() %}\n"
-        "    strt_l{{layer}} = {{value}}\n"
-        "    {%- endfor -%}"
-    )
-
-    def __init__(self, ibound, top, bottom, starting_head, inactive_head=1.0e30):
-        self._check_ibound(ibound)
-        super(__class__, self).__init__()
-        self["ibound"] = ibound
-        self["top"] = top
-        self["bottom"] = bottom
-        self["starting_head"] = starting_head
-        self["inactive_head"] = inactive_head
-
-    def _check_ibound(self, ibound):
-        if not isinstance(ibound, xr.DataArray):
-            raise TypeError("ibound must be xarray.DataArray")
-        dims = ibound.dims
-        if not (dims == ("layer", "y", "x") or dims == ("z", "y", "x")):
-            raise ValueError(
-                f'ibound dimensions must be ("layer", "y", "x") or ("z", "y", "x"),'
-                f" got instead {dims}"
-            )
-
-    def _render(self, directory, nlayer, *args, **kwargs):
-        """
-        Renders part of runfile that ends up under [bas] section.
-        """
-        d = {}
-        for varname in ("ibound", "starting_head"):
-            d[varname] = self._compose_values_layer(varname, directory, nlayer)
-        d["inactive_head"] = self["inactive_head"].values
-
-        return self._template.render(d)
-
-    def _compose_top(self, directory):
-        """
-        Composes paths to file, or gets the appropriate scalar value for
-        a top of model domain.
-
-        Parameters
-        ----------
-        directory : str
-        """
-        da = self["top"]
-        if "x" in da.coords and "y" in da.coords:
-            if not len(da.shape) == 2:
-                raise ValueError("Top should either be 2d or a scalar value")
-            d = {}
-            d["name"] = "top"
-            d["directory"] = directory
-            d["extension"] = ".idf"
-            value = self._compose(d)
-        else:
-            if not da.shape == ():
-                raise ValueError("Top should either be 2d or a scalar value")
-            value = float(da)
-        return value
-
-    @staticmethod
-    def _cellsizes(dx):
-        ncell = dx.size
-        index_ends = np.argwhere(np.diff(dx) != 0.0) + 1
-        index_ends = np.append(index_ends, ncell)
-        index_starts = np.insert(index_ends[:-1], 0, 0) + 1
-
-        d = {}
-        for s, e in zip(index_starts, index_ends):
-            value = abs(float(dx[s - 1]))
-            if s == e:
-                d[f"{s}"] = value
-            else:
-                d[f"{s}:{e}"] = value
-        return d
-
-    def _render_dis(self, directory, nlayer):
-        """
-        Renders part of runfile that ends up under [dis] section.
-        """
-        d = {}
-        d["top"] = self._compose_top(directory)
-        d["bottom"] = self._compose_values_layer("bottom", directory, nlayer)
-        d["nlay"], d["nrow"], d["ncol"] = self["ibound"].shape
-        # TODO: check dx > 0, dy < 0?
-        if "dx" not in self or "dy" not in self:  # assume equidistant
-            dx, _, _ = util.coord_reference(self["x"])
-            dy, _, _ = util.coord_reference(self["y"])
-        else:
-            dx = self.coords["dx"]
-            dy = self.coords["dy"]
-
-        if isinstance(dy, (float, int)) or dy.shape in ((), (1,)):
-            d["dy"] = {"?": abs(float(dy))}
-        else:
-            d["dy"] = self._cellsizes(dy)
-
-        if isinstance(dx, (float, int)) or dx.shape in ((), (1,)):
-            d["dx"] = {"?": abs(float(dx))}
-        else:
-            d["dx"] = self._cellsizes(dx)
-
-        # Non-time dependent part of dis
-        # Can be inferred from ibound
-        _dis_template = jinja2.Template(
-            "[dis]\n"
-            "    nlay = {{nlay}}\n"
-            "    nrow = {{nrow}}\n"
-            "    ncol = {{ncol}}\n"
-            "    {%- for row, value in dy.items() %}\n"
-            "    delc_r{{row}} = {{value}}\n"
-            "    {%- endfor %}\n"
-            "    {%- for col, value in dx.items() %}\n"
-            "    delr_c{{col}} = {{value}}\n"
-            "    {%- endfor %}\n"
-            "    top = {{top}}\n"
-            "    {%- for layer, value in bottom.items() %}\n"
-            "    botm_l{{layer}} = {{value}}\n"
-            "    {%- endfor %}\n"
-            "    laycbd_l? = 0"
-        )
-
-        return _dis_template.render(d)
-
-    def thickness(self):
-        """
-        Computes layer thickness from top and bottom data.
-
-        Returns
-        -------
-        thickness : xr.DataArray
-        """
-        th = xr.concat(
-            [
-                self["top"] - self["bottom"].sel(layer=1),
-                -1.0 * self["bottom"].diff("layer"),
-            ],
-            dim="layer",
-        )
-        return th
-
-    def _pkgcheck(self, ibound=None):
-        if (self["top"] < self["bottom"]).any():
-            raise ValueError(f"top should be larger than bottom in {self}")
-
-        active_cells = self["ibound"] != 0
-        if (active_cells & np.isnan(self["starting_head"])).any():
-            raise ValueError(
-                f"Active cells in ibound may not have a nan value in starting_head in {self}"
-            )
-
-        _, nlabels = scipy.ndimage.label(active_cells.values)
-        if nlabels > 1:
-            raise ValueError(
-                f"{nlabels} disconnected model domain detected in the ibound in {self}"
-            )

From ce5ac5eccd6a66de2ec4a51b408e159490fda972 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 14 Jun 2022 16:39:50 +0200
Subject: [PATCH 29/35] Ensure the applying the factor doesn't cast to a
 different type

---
 imod/mf6/read_input/grid_data.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/imod/mf6/read_input/grid_data.py b/imod/mf6/read_input/grid_data.py
index bee18fdba..8532f42d0 100644
--- a/imod/mf6/read_input/grid_data.py
+++ b/imod/mf6/read_input/grid_data.py
@@ -172,7 +172,8 @@ def read_array(
         )
 
     if factor is not None:
-        array = array * factor
+        # Cast the factor as well to make sure type is preserved.
+        array = array * dtype(factor)
 
     return array
 

From aa439d224f636896694db3713e239d02eadf5879 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 14 Jun 2022 16:39:58 +0200
Subject: [PATCH 30/35] Space in warning

---
 imod/mf6/read_input/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/imod/mf6/read_input/__init__.py b/imod/mf6/read_input/__init__.py
index f788c2e45..b7bcb5116 100644
--- a/imod/mf6/read_input/__init__.py
+++ b/imod/mf6/read_input/__init__.py
@@ -179,7 +179,7 @@ def tdis_time(tdis: Dict[str, Any]) -> np.ndarray:
     else:
         possibilities = ", ".join(list(TIME_UNITS.keys())[1:])
         warnings.warn(
-            "Cannot convert time coordinate to datetime."
+            "Cannot convert time coordinate to datetime. "
             "Falling back to (unitless) floating point time coordinate. \n"
             f"time_units should be one of: {possibilities}.\n"
             "start_date_time should be according to ISO 8601."

From 4abab8e03533dae5dac2e59a236779b6018ef0f6 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 14 Jun 2022 16:40:26 +0200
Subject: [PATCH 31/35] Add angledegx for connection normal in disu from dis
 method

---
 imod/mf6/disu.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
index a58ecfdda..f905f452a 100644
--- a/imod/mf6/disu.py
+++ b/imod/mf6/disu.py
@@ -103,6 +103,7 @@ def __init__(
         ihc,
         cl12,
         hwva,
+        angledegx,
         idomain=None,
     ):
         super().__init__(locals())
@@ -116,6 +117,7 @@ def __init__(
         self.dataset["ihc"] = ihc
         self.dataset["cl12"] = cl12
         self.dataset["hwva"] = hwva
+        self.dataset["angledegx"] = angledegx
         if idomain is not None:
             self.dataset["idomain"] = idomain
 
@@ -251,6 +253,7 @@ def from_dis(
         ihc = is_horizontal.astype(np.int32)
         cl12 = np.zeros_like(ihc, dtype=np.float64)
         hwva = np.zeros_like(ihc, dtype=np.float64)
+        angledegx = np.zeros_like(ihc, dtype=np.float64)
         # Fill.
         cl12[is_x] = 0.5 * dxx[index_x]  # cell center to vertical face
         cl12[is_y] = 0.5 * dyy[index_y]  # cell center to vertical face
@@ -258,6 +261,7 @@ def from_dis(
         hwva[is_x] = dyy[index_x]  # width
         hwva[is_y] = dxx[index_y]  # width
         hwva[is_vertical] = area[index_v]  # area
+        angledegx[is_y] = 90.0  # angle between connection normal and x-axis.
 
         # Set "node" and "nja" as the dimension in accordance with MODFLOW6.
         # Should probably be updated if we could move to UGRID3D...
@@ -276,6 +280,7 @@ def from_dis(
                 ihc=xr.DataArray(ihc, dims=["nja"]),
                 cl12=xr.DataArray(cl12, dims=["nja"]),
                 hwva=xr.DataArray(hwva, dims=["nja"]),
+                angledegx=xr.DataArray(angledegx, dims=["nja"]),
             )
             cell_ids = np.cumsum(active) - 1
             cell_ids[~active] = -1
@@ -292,5 +297,6 @@ def from_dis(
                 ihc=xr.DataArray(ihc, dims=["nja"]),
                 cl12=xr.DataArray(cl12, dims=["nja"]),
                 hwva=xr.DataArray(hwva, dims=["nja"]),
+                angledegx=xr.DataArray(angledegx, dims=["nja"]),
                 idomain=xr.DataArray(active.astype(np.int32), dims=["node"]),
             )

From 9b29cb1e52485c7ff140b0d0492561b7bad8bf2b Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 15 Jun 2022 15:20:49 +0200
Subject: [PATCH 32/35] Include angldegx correctly in disu input

---
 imod/mf6/disu.py               | 14 ++++++++------
 imod/templates/mf6/gwf-disu.j2 |  5 +++--
 2 files changed, 11 insertions(+), 8 deletions(-)

diff --git a/imod/mf6/disu.py b/imod/mf6/disu.py
index f905f452a..fe7fc3ac8 100644
--- a/imod/mf6/disu.py
+++ b/imod/mf6/disu.py
@@ -74,6 +74,7 @@ class LowLevelUnstructuredDiscretization(Package):
     ihc: xr.DataArray of integers
     cl12: xr.DataArray of floats
     hwva: xr.DataArray of floats
+    angldegx: xr.DataArray of floats
     """
 
     _pkg_id = "disu"
@@ -87,6 +88,7 @@ class LowLevelUnstructuredDiscretization(Package):
         "cl12": np.float64,
         "hwva": np.float64,
         "idomain": np.int32,
+        "angldegx": np.float64,
     }
     _keyword_map = {}
     _template = Package._initialize_template(_pkg_id)
@@ -103,7 +105,7 @@ def __init__(
         ihc,
         cl12,
         hwva,
-        angledegx,
+        angldegx,
         idomain=None,
     ):
         super().__init__(locals())
@@ -117,7 +119,7 @@ def __init__(
         self.dataset["ihc"] = ihc
         self.dataset["cl12"] = cl12
         self.dataset["hwva"] = hwva
-        self.dataset["angledegx"] = angledegx
+        self.dataset["angldegx"] = angldegx
         if idomain is not None:
             self.dataset["idomain"] = idomain
 
@@ -253,7 +255,7 @@ def from_dis(
         ihc = is_horizontal.astype(np.int32)
         cl12 = np.zeros_like(ihc, dtype=np.float64)
         hwva = np.zeros_like(ihc, dtype=np.float64)
-        angledegx = np.zeros_like(ihc, dtype=np.float64)
+        angldegx = np.zeros_like(ihc, dtype=np.float64)
         # Fill.
         cl12[is_x] = 0.5 * dxx[index_x]  # cell center to vertical face
         cl12[is_y] = 0.5 * dyy[index_y]  # cell center to vertical face
@@ -261,7 +263,7 @@ def from_dis(
         hwva[is_x] = dyy[index_x]  # width
         hwva[is_y] = dxx[index_y]  # width
         hwva[is_vertical] = area[index_v]  # area
-        angledegx[is_y] = 90.0  # angle between connection normal and x-axis.
+        angldegx[is_y] = 90.0  # angle between connection normal and x-axis.
 
         # Set "node" and "nja" as the dimension in accordance with MODFLOW6.
         # Should probably be updated if we could move to UGRID3D...
@@ -280,7 +282,7 @@ def from_dis(
                 ihc=xr.DataArray(ihc, dims=["nja"]),
                 cl12=xr.DataArray(cl12, dims=["nja"]),
                 hwva=xr.DataArray(hwva, dims=["nja"]),
-                angledegx=xr.DataArray(angledegx, dims=["nja"]),
+                angldegx=xr.DataArray(angldegx, dims=["nja"]),
             )
             cell_ids = np.cumsum(active) - 1
             cell_ids[~active] = -1
@@ -297,6 +299,6 @@ def from_dis(
                 ihc=xr.DataArray(ihc, dims=["nja"]),
                 cl12=xr.DataArray(cl12, dims=["nja"]),
                 hwva=xr.DataArray(hwva, dims=["nja"]),
-                angledegx=xr.DataArray(angledegx, dims=["nja"]),
+                angldegx=xr.DataArray(angldegx, dims=["nja"]),
                 idomain=xr.DataArray(active.astype(np.int32), dims=["node"]),
             )
diff --git a/imod/templates/mf6/gwf-disu.j2 b/imod/templates/mf6/gwf-disu.j2
index 2990e1aa9..82e8bdf66 100644
--- a/imod/templates/mf6/gwf-disu.j2
+++ b/imod/templates/mf6/gwf-disu.j2
@@ -41,6 +41,7 @@ begin connectiondata
     {{cl12}}
   hwva
     {{hwva}}
-{%- if angldegx is defined %}  angldegx
-    {angldegx}{% endif %}
+{% if angldegx is defined %}  angldegx
+    {{angldegx}}
+{% endif -%}
 end connectiondata
\ No newline at end of file

From 38cdafb96edcf26a7525cde612a0807292c03d55 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Mon, 20 Jun 2022 10:15:26 +0200
Subject: [PATCH 33/35] Use node number for DISU list based input, rather than
 to start numbering at node 0 for every subset of a grid.

---
 imod/mf6/pkgbase.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/imod/mf6/pkgbase.py b/imod/mf6/pkgbase.py
index e4467d25e..61b76d653 100644
--- a/imod/mf6/pkgbase.py
+++ b/imod/mf6/pkgbase.py
@@ -51,7 +51,7 @@ def disv_recarr(arrdict, layer, notnull):
     return recarr
 
 
-def disu_recarr(arrdict, layer, notnull):
+def disu_recarr(arrdict, node, notnull):
     index_spec = [("node", np.int32)]
     field_spec = [(key, np.float64) for key in arrdict]
     sparse_dtype = np.dtype(index_spec + field_spec)
@@ -59,7 +59,11 @@ def disu_recarr(arrdict, layer, notnull):
     nrow = notnull.sum()
     recarr = np.empty(nrow, dtype=sparse_dtype)
     # Argwhere returns an index_array with dims (N, a.ndims)
-    recarr["node"] = (np.argwhere(notnull) + 1)[:, 0]
+    index = np.argwhere(notnull)[:, 0]
+    if node is None:
+        recarr["node"] = index + 1
+    else:
+        recarr["node"] = node[index] + 1
     return recarr
 
 
@@ -598,6 +602,10 @@ def write_datafile(self, outpath, ds, binary):
         """
         if "layer" in ds:
             layer = ds["layer"].values
+        # TODO: change layer argument name to cell_id or something?
+        # Also forward the node number in case of DISU input.
+        elif "node" in ds:
+            layer = ds["node"].values
         else:
             layer = None
 

From 72b4842b4ed12851c53a82d1179baf2f44398bc4 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 23 Jun 2022 17:33:05 +0200
Subject: [PATCH 34/35] Add index_col=False to read_text_listinput

---
 imod/mf6/read_input/list_input.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/imod/mf6/read_input/list_input.py b/imod/mf6/read_input/list_input.py
index 7222ddbeb..090ddd351 100644
--- a/imod/mf6/read_input/list_input.py
+++ b/imod/mf6/read_input/list_input.py
@@ -53,6 +53,7 @@ def read_text_listinput(
     df = pd.read_csv(
         path,
         header=None,
+        index_col=False,
         dtype=d,
         names=d.keys(),
         delim_whitespace=True,

From fe172a5eac48ed30c8bcb06fd929c1fd00b5db11 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 17 Aug 2022 09:53:10 +0200
Subject: [PATCH 35/35] Take 1-based indexing into account for MODFLOW periods

---
 imod/mf6/sto.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/imod/mf6/sto.py b/imod/mf6/sto.py
index ec957afc8..583646bce 100644
--- a/imod/mf6/sto.py
+++ b/imod/mf6/sto.py
@@ -65,9 +65,10 @@ def open(
         for field, data in griddata.items():
             content[field] = xr.DataArray(data, coords, dims)
         periods = content.pop("periods")
+        time_index = np.fromiter(periods.keys(), dtype=int) - 1
         content["transient"] = xr.DataArray(
             list(periods.values()),
-            coords={"time": globaltimes[list(periods.keys())]},
+            coords={"time": globaltimes[time_index]},
             dims=("time",),
         )