From c4c3d2ad364783db7f099fe4d0afbd33119bc170 Mon Sep 17 00:00:00 2001
From: Eric Morway <emorway@usgs.gov>
Date: Wed, 22 May 2024 16:16:09 -0700
Subject: [PATCH] stop writing blank entries to output, patch up some autotests
 as a result

---
 autotest/test_gwf_mvr01.py            |  57 ++++--------
 autotest/test_gwt_uztmvt2x2.py        |  13 ++-
 src/Model/GroundWaterFlow/gwf-sfr.f90 | 120 ++++++++++++++------------
 3 files changed, 88 insertions(+), 102 deletions(-)

diff --git a/autotest/test_gwf_mvr01.py b/autotest/test_gwf_mvr01.py
index 0913dbc40c2..78f84b18988 100644
--- a/autotest/test_gwf_mvr01.py
+++ b/autotest/test_gwf_mvr01.py
@@ -379,48 +379,29 @@ def check_output(idx, test):
     times = bobj.get_times()
     records = bobj.get_data(totim=times[-1])
     adt = [("node", "<i4"), ("node2", "<i4"), ("q", "<f8")]
-    assert len(records) == 25
-    assert records[0].shape == (0,)
+    assert len(records) == 5
 
-    assert records[1].shape == (2,)
+    assert records[0].shape == (2,)
     a = np.array([(1, 2, -0.0), (2, 2, -0.0)], dtype=adt)
-    assert np.array_equal(records[1], a)
+    assert np.array_equal(records[0], a)
 
-    assert records[2].shape == (2,)
+    assert records[1].shape == (2,)
     a = np.array([(1, 1, -0.0), (2, 1, -0.0)], dtype=adt)
-    assert np.array_equal(records[2], a)
+    assert np.array_equal(records[1], a)
 
-    assert records[3].shape == (2,)
+    assert records[2].shape == (2,)
     a = np.array([(1, 3, -0.00545875), (2, 3, -0.00468419)], dtype=adt)
-    assert np.allclose(records[3]["node"], a["node"])
-    assert np.allclose(records[3]["node2"], a["node2"])
-    assert np.allclose(records[3]["q"], a["q"], atol=0.001), "{}\n{}".format(
-        records[3]["q"], a["q"]
+    assert np.allclose(records[2]["node"], a["node"])
+    assert np.allclose(records[2]["node2"], a["node2"])
+    assert np.allclose(records[2]["q"], a["q"], atol=0.001), "{}\n{}".format(
+        records[2]["q"], a["q"]
     )
 
-    assert records[4].shape == (0,)
-    assert records[5].shape == (0,)
-    assert records[6].shape == (0,)
-    assert records[7].shape == (0,)
-    assert records[8].shape == (3,)
+    assert records[3].shape == (3,)
     a = np.array([(1, 1, -0.0), (1, 2, -0.0005), (1, 3, -0.0)], dtype=adt)
-    assert np.array_equal(records[8], a)
-
-    assert records[9].shape == (0,)
-    assert records[10].shape == (0,)
-    assert records[11].shape == (0,)
-    assert records[12].shape == (0,)
-    assert records[13].shape == (0,)
-    assert records[14].shape == (0,)
-    assert records[15].shape == (0,)
-    assert records[16].shape == (0,)
-    assert records[17].shape == (0,)
-    assert records[18].shape == (0,)
-    assert records[19].shape == (0,)
-    assert records[20].shape == (0,)
-    assert records[21].shape == (0,)
-    assert records[22].shape == (0,)
-    assert records[23].shape == (9,)
+    assert np.array_equal(records[3], a)
+
+    assert records[4].shape == (9,)
     a = np.array(
         [
             (1, 1, -1.0e-04),
@@ -435,14 +416,12 @@ def check_output(idx, test):
         ],
         dtype=adt,
     )
-    assert np.allclose(records[23]["node"], a["node"])
-    assert np.allclose(records[23]["node2"], a["node2"])
-    assert np.allclose(records[23]["q"], a["q"], atol=0.001), "{}\n{}".format(
-        records[23]["q"], a["q"]
+    assert np.allclose(records[4]["node"], a["node"])
+    assert np.allclose(records[4]["node2"], a["node2"])
+    assert np.allclose(records[4]["q"], a["q"], atol=0.001), "{}\n{}".format(
+        records[4]["q"], a["q"]
     )
 
-    assert records[24].shape == (0,)
-
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))
 def test_mf6model(idx, name, function_tmpdir, targets):
diff --git a/autotest/test_gwt_uztmvt2x2.py b/autotest/test_gwt_uztmvt2x2.py
index c3d12ed8f78..f9e90da5ced 100644
--- a/autotest/test_gwt_uztmvt2x2.py
+++ b/autotest/test_gwt_uztmvt2x2.py
@@ -622,14 +622,11 @@ def check_output(idx, test):
 
     msg2 = "Mass received by SFR ('FROM-MVR') not as expected"
     for x in np.arange(len(fromMvrDat)):
-        for y in np.arange(len(fromMvrDat[x + 1][0])):
-            if fromMvrDat[x + 1][0][y][-1] == concCell[y]:
-                continue
-            else:
-                for z in np.arange(len(mvtdat[x + 1][y])):
-                    assert np.isclose(
-                        mvtdat[x + 1][y][z][-1], (x + 1.0) * concCell[z]
-                    ), msg2
+        for y in np.arange(len(fromMvrDat[x + 1])):
+            for z in np.arange(len(mvtdat[x + 1][y])):
+                assert np.isclose(
+                    mvtdat[x + 1][y][z][-1], (x + 1.0) * concCell[z]
+                ), msg2
 
 
 # - No need to change any code below
diff --git a/src/Model/GroundWaterFlow/gwf-sfr.f90 b/src/Model/GroundWaterFlow/gwf-sfr.f90
index b8333475b26..e5b1d8733b7 100644
--- a/src/Model/GroundWaterFlow/gwf-sfr.f90
+++ b/src/Model/GroundWaterFlow/gwf-sfr.f90
@@ -2442,8 +2442,7 @@ end subroutine sfr_cq
 
   !> @ brief Output package flow terms.
     !!
-    !!  Output SFR package flow terms.
-    !!
+    !! Write flows to binary file and/or print flows to budget
   !<
   subroutine sfr_ot_package_flows(this, icbcfl, ibudfl)
     ! -- modules
@@ -2451,7 +2450,7 @@ subroutine sfr_ot_package_flows(this, icbcfl, ibudfl)
     ! -- dummy variables
     class(SfrType) :: this !< SfrType object
     integer(I4B), intent(in) :: icbcfl !< flag and unit number for cell-by-cell output
-    integer(I4B), intent(in) :: ibudfl !< flag indication if cell-by-cell data should be saved
+    integer(I4B), intent(in) :: ibudfl !< flag indicating if cell-by-cell data should be saved
     ! -- local variables
     integer(I4B) :: ibinun
     character(len=20), dimension(:), allocatable :: cellidstr
@@ -5233,6 +5232,11 @@ subroutine sfr_setup_budobj(this)
     !    so they can be written to the binary budget files, but these internal
     !    flows are not included as part of the budget table.
     nbudterm = 8
+    !
+    ! -- GWF models with a single reach should not have an entry
+    if (this%nconn == 0) nbudterm = nbudterm - 1
+    !
+    ! -- Account for mover connection and aux vars in budget object
     if (this%imover == 1) nbudterm = nbudterm + 2
     if (this%naux > 0) nbudterm = nbudterm + 1
     !
@@ -5243,29 +5247,33 @@ subroutine sfr_setup_budobj(this)
     idx = 0
     !
     ! -- Go through and set up each budget term
-    text = '    FLOW-JA-FACE'
-    idx = idx + 1
-    maxlist = this%nconn
-    naux = 1
-    auxtxt(1) = '       FLOW-AREA'
-    call this%budobj%budterm(idx)%initialize(text, &
-                                             this%name_model, &
-                                             this%packName, &
-                                             this%name_model, &
-                                             this%packName, &
-                                             maxlist, .false., .false., &
-                                             naux, auxtxt)
     !
-    ! -- store connectivity
-    call this%budobj%budterm(idx)%reset(this%nconn)
-    q = DZERO
-    do n = 1, this%maxbound
-      n1 = n
-      do i = this%ia(n) + 1, this%ia(n + 1) - 1
-        n2 = this%ja(i)
-        call this%budobj%budterm(idx)%update_term(n1, n2, q)
+    ! -- GWF models with a single reach do not have flow-ja-face connections
+    if (this%nconn /= 0) then
+      text = '    FLOW-JA-FACE'
+      idx = idx + 1
+      maxlist = this%nconn
+      naux = 1
+      auxtxt(1) = '       FLOW-AREA'
+      call this%budobj%budterm(idx)%initialize(text, &
+                                               this%name_model, &
+                                               this%packName, &
+                                               this%name_model, &
+                                               this%packName, &
+                                               maxlist, .false., .false., &
+                                               naux, auxtxt)
+      !
+      ! -- store connectivity
+      call this%budobj%budterm(idx)%reset(this%nconn)
+      q = DZERO
+      do n = 1, this%maxbound
+        n1 = n
+        do i = this%ia(n) + 1, this%ia(n + 1) - 1
+          n2 = this%ja(i)
+          call this%budobj%budterm(idx)%update_term(n1, n2, q)
+        end do
       end do
-    end do
+    end if
     !
     ! --
     text = '             GWF'
@@ -5456,40 +5464,42 @@ subroutine sfr_fill_budobj(this)
     idx = 0
     !
     ! -- FLOW JA FACE
-    idx = idx + 1
-    call this%budobj%budterm(idx)%reset(this%nconn)
-    do n = 1, this%maxbound
-      n1 = n
-      q = DZERO
-      ca = DZERO
-      do i = this%ia(n) + 1, this%ia(n + 1) - 1
-        n2 = this%ja(i)
-        if (this%iboundpak(n) /= 0) then
-          ! flow to downstream reaches
-          if (this%idir(i) < 0) then
-            qt = this%dsflow(n)
-            q = -this%qconn(i)
-            ! flow from upstream reaches
+    if (this%nconn /= 0) then
+      idx = idx + 1
+      call this%budobj%budterm(idx)%reset(this%nconn)
+      do n = 1, this%maxbound
+        n1 = n
+        q = DZERO
+        ca = DZERO
+        do i = this%ia(n) + 1, this%ia(n + 1) - 1
+          n2 = this%ja(i)
+          if (this%iboundpak(n) /= 0) then
+            ! flow to downstream reaches
+            if (this%idir(i) < 0) then
+              qt = this%dsflow(n)
+              q = -this%qconn(i)
+              ! flow from upstream reaches
+            else
+              qt = this%usflow(n)
+              do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
+                if (this%idir(ii) > 0) cycle
+                if (this%ja(ii) /= n) cycle
+                q = this%qconn(ii)
+                exit
+              end do
+            end if
+            ! calculate flow area
+            call this%sfr_calc_reach_depth(n, qt, d)
+            ca = this%calc_area_wet(n, d)
           else
-            qt = this%usflow(n)
-            do ii = this%ia(n2) + 1, this%ia(n2 + 1) - 1
-              if (this%idir(ii) > 0) cycle
-              if (this%ja(ii) /= n) cycle
-              q = this%qconn(ii)
-              exit
-            end do
+            q = DZERO
+            ca = DZERO
           end if
-          ! calculate flow area
-          call this%sfr_calc_reach_depth(n, qt, d)
-          ca = this%calc_area_wet(n, d)
-        else
-          q = DZERO
-          ca = DZERO
-        end if
-        this%qauxcbc(1) = ca
-        call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc)
+          this%qauxcbc(1) = ca
+          call this%budobj%budterm(idx)%update_term(n1, n2, q, this%qauxcbc)
+        end do
       end do
-    end do
+    end if
     !
     ! -- GWF (LEAKAGE)
     idx = idx + 1