implement fix for lake status set to inactive

autotest/test_gwf_lak_status.py

@@ -198,7 +198,7 @@ def check_budget_file(idx, test):
  print (f"Checking contents of {fpth.name}")
  bobj = flopy.utils.CellBudgetFile(fpth, precision="double")
  times = bobj.get_times()
- # bobj.list_unique_records()
+ print (bobj.list_unique_records())
 
  # check to make sure GWF is zero when lake is inactive
  # and non-zero otherwise
@@ -211,11 +211,23 @@ def check_budget_file(idx, test):
  else:
  assert r["q"] != 0.
 
+ # check all the other terms to make sure they are zero for 
+ # the second period
+ lake_terms = [
+ "rainfall", "evaporation", "ext-inflow", "withdrawal", 
+ "ext-outflow", "storage", "constant"
+ ]
+ for lake_term in lake_terms:
+ print (f" Checking lake term {lake_term} for period 2.")
+ record = bobj.get_data(text=lake_term, totim=times[1])[0]
+ q = record["q"][0]
+ msg = f"Lake term for stress period 2 is not zero ({lake_term}={q})"
+ assert np.allclose(q, 0.), msg
 
 def check_stage_file(idx, test):
- # Check to make sure the stage is equal to dnodata
- # for stress period 2 and not dnodata otherwise
- dnodata = 1.e30
+ # Check to make sure the stage is equal to dhnoflo
+ # for stress period 2 and not dhnoflo otherwise
+ dhnoflo = 1.e30
  name = cases[idx]
  fpth = test.workspace / f"{name}.lak.stage"
  print (f"Checking contents of {fpth.name}")
@@ -224,35 +236,50 @@ def check_stage_file(idx, test):
  for kper in range(3):
  print (f" Checking binary stage for stress period {kper + 1}")
  stage = bobj.get_data(totim=times[kper]).flatten()
+ print (stage)
  if kper == 1:
- assert stage[0] == dnodata
+ assert stage[0] == dhnoflo
  else:
- assert stage[0] != dnodata
+ assert stage[0] != dhnoflo
 
 
 def check_head_file(idx, test):
  # Check to make sure the simulated head for the first and 
  # last period is the same and different from the second period.
- dnodata = 1.e30
  name = cases[idx]
  fpth = test.workspace / f"{name}.hds"
  print (f"Checking contents of {fpth.name}")
  bobj = flopy.utils.HeadFile(fpth, text="head", precision="double")
  times = bobj.get_times()
 
  head0 = bobj.get_data(totim=times[0]).flatten()
- head1 = bobj.get_data(totim=times[0]).flatten()
- head2 = bobj.get_data(totim=times[0]).flatten()
+ head1 = bobj.get_data(totim=times[1]).flatten()
+ head2 = bobj.get_data(totim=times[2]).flatten()
 
  assert np.allclose(head0, head2), (
  "Simulated heads for period 1 and 3 should be the same"
  )
  assert np.all(np.less_equal(head1, 100.)), (
- "Simulated heads for period 2 should be less than or equal to 100."
+ f"Simulated heads for period 2 should be less than or equal to 100. {head1}"
  )
 
 
+def check_lake_obs(idx, test):
+ # Check to make sure the simulated head for the first and 
+ # last period is the same and different from the second period.
+ dnodata = 3.e30
+ name = cases[idx]
+ fpth = test.workspace / f"{name}.lak.obs.csv"
+ print (f"Checking contents of {fpth.name}")
+ obs = np.genfromtxt(fpth, names=True, delimiter=",")
+ stage = obs["LAKESTAGE"].tolist()
+ print (stage)
+ assert stage[0] == stage[2], "Period 1 and period 3 stages should be equal."
+ assert stage[1] == dnodata, "Period 2 stage should equal dnodata"
+
+
 def check_output(idx, test):
+ check_lake_obs(idx, test)
  check_head_file(idx, test)
  check_stage_file(idx, test)
  check_budget_file(idx, test)

src/Model/GroundWaterFlow/gwf-lak.f90

@@ -238,7 +238,6 @@ module LakModule
  procedure, private :: lak_check_valid
  procedure, private :: lak_set_stressperiod
  procedure, private :: lak_set_attribute_error
- procedure, private :: lak_cfupdate
  procedure, private :: lak_bound_update
  procedure, private :: lak_calculate_sarea
  procedure, private :: lak_calculate_warea
@@ -3610,9 +3609,6 @@ subroutine lak_cf(this)
  integer(I4B) :: igwfnode
  real(DP) :: hlak, blak
  !
- ! -- Calculate lak conductance and update package RHS and HCOF
- !call this%lak_cfupdate()
- !
  ! -- save groundwater seepage for lake solution
  do n = 1, this%nlakes
  this%seep0(n) = this%seep(n)
@@ -3711,6 +3707,7 @@ subroutine lak_fc(this, rhs, ia, idxglo, matrix_sln)
  !
  ! -- add terms to the gwf matrix
  do n = 1, this%nlakes
+ if (this%iboundpak(n) == 0) cycle
  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
  igwfnode = this%cellid(j)
  if (this%ibound(igwfnode) < 1) cycle
@@ -4170,15 +4167,16 @@ subroutine lak_cq(this, x, flowja, iadv)
  !
  ! -- gwf flow and constant flow to lake
  do n = 1, this%nlakes
- if (this%iboundpak(n) == 0) cycle
  rrate = DZERO
  do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
  ! simvals is from aquifer perspective, and so it is positive
  ! for flow into the aquifer. Need to switch sign for lake
  ! perspective.
  rrate = -this%simvals(j)
  this%qleak(j) = rrate
- call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
+ if (this%iboundpak(n) /= 0) then
+ call this%lak_accumulate_chterm(n, rrate, chratin, chratout)
+ end if
  end do
  end do
  !
@@ -4282,10 +4280,17 @@ subroutine lak_ot_dv(this, idvsave, idvprint)
  !
  ! -- write data
  do n = 1, this%nlakes
- stage = this%xnewpak(n)
- call this%lak_calculate_sarea(n, stage, sa)
- call this%lak_calculate_warea(n, stage, wa)
- call this%lak_calculate_vol(n, stage, v)
+ if (this%iboundpak(n) == 0) then
+ stage = DHNOFLO
+ sa = DHNOFLO
+ wa = DHNOFLO
+ v = DHNOFLO
+ else
+ stage = this%xnewpak(n)
+ call this%lak_calculate_sarea(n, stage, sa)
+ call this%lak_calculate_warea(n, stage, wa)
+ call this%lak_calculate_vol(n, stage, v)
+ end if
  if (this%inamedbound == 1) then
  call this%stagetab%add_term(this%lakename(n))
  end if
@@ -5070,51 +5075,6 @@ subroutine lak_accumulate_chterm(this, ilak, rrate, chratin, chratout)
  return
  end subroutine lak_accumulate_chterm
 
- !> @brief Update LAK satcond and package rhs and hcof
- !<
- subroutine lak_cfupdate(this)
- ! -- dummy
- class(LakType), intent(inout) :: this
- ! -- local
- integer(I4B) :: j, n, node
- real(DP) :: hlak, head, clak, blak
- !
- ! -- Return if no lak lakes
- if (this%nbound .eq. 0) return
- !
- ! -- Calculate hcof and rhs for each lak entry
- do n = 1, this%nlakes
- hlak = this%xnewpak(n)
- do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
- node = this%cellid(j)
- head = this%xnew(node)
-
- this%hcof(j) = DZERO
- this%rhs(j) = DZERO
- !
- ! -- set bound, hcof, and rhs components
- call this%lak_calculate_conn_conductance(n, j, hlak, head, clak)
- this%simcond(j) = clak
-
- this%bound(2, j) = clak
-
- blak = this%bound(3, j)
-
- this%hcof(j) = -clak
- !
- ! -- fill rhs
- if (hlak < blak) then
- this%rhs(j) = -clak * blak
- else
- this%rhs(j) = -clak * hlak
- end if
- end do
- end do
- !
- ! -- Return
- return
- end subroutine lak_cfupdate
-
  !> @brief Store the lake head and connection conductance in the bound array
  !<
  subroutine lak_bound_update(this)
@@ -5281,9 +5241,13 @@ subroutine lak_solve(this, update)
  this%flwiter(n) = this%flwin(n)
  this%flwiter1(n) = this%flwin(n)
  if (this%gwfiss /= 0) then
- this%flwiter(n) = DEP20 !1.D+10
- this%flwiter1(n) = DEP20 !1.D+10
+ this%flwiter(n) = DEP20
+ this%flwiter1(n) = DEP20
  end if
+ do j = this%idxlakeconn(n), this%idxlakeconn(n + 1) - 1
+ this%hcof(j) = DZERO
+ this%rhs(j) = DZERO
+ end do
  end do
  !
  estseep: do i = 1, 2
@@ -5335,6 +5299,11 @@ subroutine lak_solve(this, update)
  end do estseep
  !
  laklevel: do n = 1, this%nlakes
+ ! -- skip inactive lakes
+ if (this%iboundpak(n) == 0) then
+ this%ncncvr(n) = 1
+ cycle laklevel
+ end if
  ibflg = 0
  hlak = this%xnewpak(n)
  if (iter < maxiter) then