write structured grid arrays 1 row per line

flopy4/mf6/codec/writer/templates/macros.jinja

@@ -44,13 +44,13 @@
 
 {{ inset }}INTERNAL
 {% for chunk in layer|array2chunks -%}
-{{ (2 * inset) ~ chunk|array2string(context.float_precision) }}
+{{ (2 * inset) ~ chunk|array2string(context.float_precision)|indent(2 * inset|length) }}
 {%- endfor %}
 {%- endfor %}
 {% elif how == "internal" %}
 {{ inset }}INTERNAL
 {% for chunk in value|array2chunks -%}
-{{ (2 * inset) ~ chunk|array2string(context.float_precision) }}
+{{ (2 * inset) ~ chunk|array2string(context.float_precision)|indent(2 * inset|length) }}
 {%- endfor %}
 {% elif how == "external" %}
 OPEN/CLOSE {{ value|quote_if_needed }}

flopy4/mf6/converter/egress/unstructure.py

@@ -298,25 +298,35 @@ def _unstructure_package(value: Package) -> dict[str, Any]:
             if meta.get("index"):
                 field_value = _to_file_index(field_value)
             if meta["shape"]:
-                # reshape layered array to (nlay, ncpl) with named
-                # dims to signal the writer to use layered format
-                if (
-                    meta.get("layered")
-                    and hasattr(field_value, "reshape")
-                    and not isinstance(field_value, xr.DataArray)
-                ):
+                # shape a per-cell array like the grid: a layered one by
+                # layer, with an nlay dim to signal the writer to use
+                # layered format, and on a structured grid each layer
+                # (nrow, ncol), since MF6 reads such a layer a row per line
+                if hasattr(field_value, "reshape") and not isinstance(field_value, xr.DataArray):
                     _get_dims = getattr(value, "get_dims", None)
                     _dims_d = _get_dims() if _get_dims else {}
                     _nlay = _dims_d.get("nlay", 0)
                     _ncpl = _dims_d.get("ncpl", 0)
-                    if _nlay > 1 and _ncpl > 0 and field_value.size == _nlay * _ncpl:
-                        blocks[block_name][f.name] = _mark_netcdf(
-                            xr.DataArray(
-                                field_value.reshape(_nlay, _ncpl),
-                                dims=("nlay", "ncpl"),
-                            ),
-                            meta,
+                    _layer = (
+                        (_dims_d["nrow"], _dims_d["ncol"])
+                        if "nrow" in _dims_d and "ncol" in _dims_d
+                        else (_ncpl,)
+                    )
+                    _da = None
+                    if (
+                        meta.get("layered")
+                        and _nlay > 1
+                        and _ncpl > 0
+                        and field_value.size == _nlay * _ncpl
+                    ):
+                        _da = xr.DataArray(
+                            field_value.reshape(_nlay, *_layer),
+                            dims=("nlay", *(("nrow", "ncol") if len(_layer) > 1 else ("ncpl",))),
                         )
+                    elif len(_layer) > 1 and _ncpl > 0 and field_value.size == _ncpl:
+                        _da = xr.DataArray(field_value.reshape(_layer), dims=("nrow", "ncol"))
+                    if _da is not None:
+                        blocks[block_name][f.name] = _mark_netcdf(_da, meta)
                         continue
                 blocks[block_name][f.name] = _mark_netcdf(_wrap_array(field_value), meta)
 

test/mf6/test_mf6_integration.py

@@ -2919,3 +2919,45 @@ def test_rcha_tas_reference_netcdf(function_tmpdir):
     # NetCDF input needs the extended build
     if os.getenv("MF6_EXTENDED"):
         sim.run()
+
+
+def test_structured_grid_arrays_written_by_row(function_tmpdir):
+    """MF6 reads a structured grid's layer a row per line, so a varying
+    griddata array is written that way: TOP, a LAYERED BOTM, and a
+    non-layered K spanning both layers."""
+    sim = Simulation(
+        name="rows",
+        workspace=function_tmpdir,
+        tdis=Time(perlen=[1.0], nstp=[1]),
+        solutiongroup=[
+            Ims(
+                models=["gwf"],
+                outer_dvclose=1e-9,
+                outer_maximum=100,
+                inner_maximum=300,
+                inner_dvclose=1e-9,
+                linear_acceleration="cg",
+            )
+        ],
+    )
+    dis = Dis(
+        nlay=2,
+        nrow=3,
+        ncol=4,
+        delr=1.0,
+        delc=1.0,
+        top=np.linspace(10.0, 11.0, 12),
+        botm=np.repeat([5.0, 0.0], 12) + np.tile(np.linspace(0, 0.3, 12), 2),
+    )
+    gwf = Gwf(parent=sim, name="gwf", dis=dis)
+    Ic(parent=gwf, strt=10.0)
+    Npf(parent=gwf, k=np.linspace(1.0, 2.0, 24), icelltype=0)
+    Chd(parent=gwf, stress_period_data={0: [[(0, 0, 0), 10.0], [(1, 2, 3), 2.0]]})
+    Oc(parent=gwf, head_file="gwf.hds", stress_period_data={0: [("SAVE", "HEAD", "LAST")]})
+    sim.write()
+
+    lines = (function_tmpdir / "gwf.dis").read_text().splitlines()
+    top = lines.index(" TOP")
+    assert lines[top + 1] == " INTERNAL"
+    assert all(len(ln.split()) == 4 for ln in lines[top + 2 : top + 5])
+    sim.run()