feat(disu-idomain): write idomain to binary grid file for disu models (…

…MODFLOW-USGS#2144)

The binary grid file written by MODFLOW 6 for DISU models did not include the IDOMAIN array. The binary grid file now includes IDOMAIN for all discretization types, including DISU.

autotest/test_gwf_disu.py

@@ -47,7 +47,7 @@ def build_models(idx, test):
     ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0)
     npf = flopy.mf6.ModflowGwfnpf(gwf)
     spd = {0: [[(0,), 1.0], [(nrow * ncol - 1), 0.0]]}
-    chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd)
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=spd)
     return sim, None
 
 
@@ -59,13 +59,18 @@ def check_output(idx, test):
     nodes = grbobj._datadict["NODES"]
     ia = grbobj._datadict["IA"]
     ja = grbobj._datadict["JA"]
+    idomain = grbobj._datadict["IDOMAIN"]
+
+    if idx == 0:
+        assert np.array_equal(idomain, np.array(27 * [1]), int)
 
     if idx == 1:
         assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7]))
         assert np.array_equal(ja[:6], np.array([1, 4, 10, 3, 6, 12]))
         assert ia[-1] == 127
         assert ia.shape[0] == 28, "ia should have size of 28"
         assert ja.shape[0] == 126, "ja should have size of 126"
+        assert np.array_equal(idomain, np.array([1, 0] + 25 * [1]), int)
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))

doc/ReleaseNotes/develop.tex

@@ -6,7 +6,7 @@
 %\underline{NEW FUNCTIONALITY}
 %\begin{itemize}
 	\item Support for adjusting time step lengths using the adaptive time stepping (ATS) capability was added to the GWT Advection (ADV) Package of the Groundwater Transport (GWT) Model in release 6.6.0.  The same functionality that was added to GWT is now available with the Groundwater Energy Transport (GWE) Model.  A description of how this functionality works and how to activate it can be found in the release notes for version 6.6.0 (Appendix A) and in the MODFLOW 6 input-output guide.
-%	\item xxx
+	\item The binary grid file written by MODFLOW 6 for DISU models did not include the IDOMAIN array.  The binary grid file now includes IDOMAIN for all discretization types, including DISU.
 %	\item xxx
 %\end{itemize}
 

doc/mf6io/framework/binaryoutput.tex

@@ -157,15 +157,16 @@ \subsubsection{DISU Grids}
 \noindent Definition 7: \texttt{`BOT DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
 \noindent Definition 8: \texttt{`IA INTEGER NDIM 1 ncells+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
 \noindent Definition 9: \texttt{`JA INTEGER NDIM 1 nja'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
-\noindent Definition 10: \texttt{`ICELLTYPE INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 10: \texttt{`IDOMAIN INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 11: \texttt{`ICELLTYPE INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
 
 \vspace{5mm}
 \noindent If vertices are provided in the DISU Package, then 5 additional definitions are included: \\
-\noindent Definition 11: \texttt{`VERTICES DOUBLE NDIM 2 2 nvert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
-\noindent Definition 12: \texttt{`CELLX DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
-\noindent Definition 13: \texttt{`CELLY DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
-\noindent Definition 14: \texttt{`IAVERT INTEGER NDIM 1 nodes+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
-\noindent Definition 15: \texttt{`JAVERT INTEGER NDIM 1 njavert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 12: \texttt{`VERTICES DOUBLE NDIM 2 2 nvert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 13: \texttt{`CELLX DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 14: \texttt{`CELLY DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 15: \texttt{`IAVERT INTEGER NDIM 1 nodes+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
+\noindent Definition 16: \texttt{`JAVERT INTEGER NDIM 1 njavert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\
 
 \vspace{5mm}
 \noindent Read \texttt{NDAT} data variables using the definitions defined above. \\
@@ -178,15 +179,16 @@ \subsubsection{DISU Grids}
 \noindent Record 7: \texttt{((BOT(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\
 \noindent Record 8: \texttt{(IA(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\
 \noindent Record 9: \texttt{(JA(J),J=1,NJA)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJA)}} \\
-\noindent Record 10: \texttt{(ICELLTYPE(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\
+\noindent Record 10: \texttt{(IDOMAIN(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\
+\noindent Record 11: \texttt{(ICELLTYPE(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\
 
 \vspace{5mm}
 \noindent If vertices are provided in the DISU Package, then 5 additional records are included: \\
-\noindent Record 11: \texttt{((VERT(J,K),J=1,2),K=1,NVERT)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(2,NVERT)}} \\
-\noindent Record 12: \texttt{(CELLX(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}}\\
-\noindent Record 13: \texttt{(CELLY(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\
-\noindent Record 14: \texttt{(IAVERT(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\
-\noindent Record 15: \texttt{(JAVERT(J),J=1,NJAVERT)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJAVERT)}} \\
+\noindent Record 12: \texttt{((VERT(J,K),J=1,2),K=1,NVERT)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(2,NVERT)}} \\
+\noindent Record 13: \texttt{(CELLX(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}}\\
+\noindent Record 14: \texttt{(CELLY(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\
+\noindent Record 15: \texttt{(IAVERT(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\
+\noindent Record 16: \texttt{(JAVERT(J),J=1,NJAVERT)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJAVERT)}} \\
 
 
 \newpage

src/Model/Discretization/Disu.f90

@@ -923,7 +923,7 @@ subroutine write_grb(this, icelltype)
        &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)"
     !
     ! -- Initialize
-    ntxt = 10
+    ntxt = 11
     if (this%nvert > 0) ntxt = ntxt + 5
     !
     ! -- Open the file
@@ -975,6 +975,9 @@ subroutine write_grb(this, icelltype)
     write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', this%con%nja
     txt(lentxt:lentxt) = new_line('a')
     write (iunit) txt
+    write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
+    txt(lentxt:lentxt) = new_line('a')
+    write (iunit) txt
     write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser
     txt(lentxt:lentxt) = new_line('a')
     write (iunit) txt
@@ -1008,6 +1011,7 @@ subroutine write_grb(this, icelltype)
     write (iunit) this%bot1d ! bot
     write (iunit) this%con%iausr ! ia
     write (iunit) this%con%jausr ! ja
+    write (iunit) this%idomain ! idomain
     write (iunit) icelltype ! icelltype
     !
     ! -- if vertices have been read then write additional data