feat(idm): update RIV package for IDM (MODFLOW-USGS#1378)

* update RIV package for IDM * rebuild makefiles * cleanup --------- Co-authored-by: mjreno <mreno@IGSAAA071L00066.gs.doi.net>
emorway-usgs · Sep 30, 2023 · 135e1a6 · 135e1a6
1 parent 25e3efa
commit 135e1a6
Show file tree

Hide file tree

Showing 12 changed files with 785 additions and 101 deletions.
diff --git a/autotest/test_gwf_libmf6_riv01.py b/autotest/test_gwf_libmf6_riv01.py
@@ -178,8 +178,12 @@ def api_func(exe, idx, model_ws=None):
     end_time = mf6.get_end_time()
 
     # get copy of (multi-dim) array with river parameters
-    riv_tag = mf6.get_var_address("BOUND", name, riv_packname)
-    new_spd = mf6.get_value(riv_tag)
+    stage_tag = mf6.get_var_address("STAGE", name, riv_packname)
+    cond_tag = mf6.get_var_address("COND", name, riv_packname)
+    rbot_tag = mf6.get_var_address("RBOT", name, riv_packname)
+    new_stage = mf6.get_value(stage_tag)
+    new_cond = mf6.get_value(cond_tag)
+    new_rbot = mf6.get_value(rbot_tag)
 
     # model time loop
     idx = 0
@@ -192,16 +196,18 @@ def api_func(exe, idx, model_ws=None):
         mf6.prepare_time_step(dt)
 
         # set the RIV data through the BMI
+        # change cond and rbot data
+        new_cond[:] = [riv_cond]
+        new_rbot[:] = [riv_bot]
+        mf6.set_value(cond_tag, new_cond)
+        mf6.set_value(rbot_tag, new_rbot)
+        # change stage data
         if current_time < 5:
-            # set columns of BOUND data (we're setting entire columns of the
-            # 2D array for convenience, setting only the value for the active
-            # stress period should work too)
-            new_spd[:] = [riv_stage, riv_cond, riv_bot]
-            mf6.set_value(riv_tag, new_spd)
+            new_stage[:] = [riv_stage]
+            mf6.set_value(stage_tag, new_stage)
         else:
-            # change only stage data
-            new_spd[:] = [riv_stage2, riv_cond, riv_bot]
-            mf6.set_value(riv_tag, new_spd)
+            new_stage[:] = [riv_stage2]
+            mf6.set_value(stage_tag, new_stage)
 
         kiter = 0
         mf6.prepare_solve()

diff --git a/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn b/doc/mf6io/mf6ivar/dfn/gwf-riv.dfn
@@ -35,6 +35,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'river'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -43,6 +44,7 @@ reader urword
 optional true
 longname print calculated flows to listing file
 description REPLACE print_flows {'{#1}': 'river'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -51,6 +53,7 @@ reader urword
 optional true
 longname save CHD flows to budget file
 description REPLACE save_flows {'{#1}': 'river'}
+mf6internal ipakcb
 
 block options
 name ts_filerecord
@@ -169,6 +172,7 @@ shape (maxbound)
 reader urword
 longname
 description
+mf6internal spd
 
 block period
 name cellid
@@ -224,6 +228,7 @@ optional true
 time_series true
 longname auxiliary variables
 description REPLACE aux {'{#1}': 'river'}
+mf6internal auxvar
 
 block period
 name boundname

diff --git a/make/makefile b/make/makefile
@@ -100,6 +100,7 @@ $(OBJDIR)/gwt1disv1idm.o \
 $(OBJDIR)/gwt1disu1idm.o \
 $(OBJDIR)/gwt1dis1idm.o \
 $(OBJDIR)/gwf3wel8idm.o \
+$(OBJDIR)/gwf3riv8idm.o \
 $(OBJDIR)/gwf3npf8idm.o \
 $(OBJDIR)/gwf3idm.o \
 $(OBJDIR)/gwf3drn8idm.o \

diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj
@@ -146,6 +146,7 @@
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3oc8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sfr8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sto8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3tvbase8.f90"/>

diff --git a/src/Model/GroundWaterFlow/gwf3.f90 b/src/Model/GroundWaterFlow/gwf3.f90
@@ -1299,7 +1299,8 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
       call drn_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
                       pakname, mempath)
     case ('RIV6')
-      call riv_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
+      call riv_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
+                      pakname, mempath)
     case ('GHB6')
       call ghb_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
     case ('RCH6')

diff --git a/src/Model/GroundWaterFlow/gwf3buy8.f90 b/src/Model/GroundWaterFlow/gwf3buy8.f90
@@ -648,6 +648,7 @@ subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, &
 ! ------------------------------------------------------------------------------
     ! -- modules
     use BndModule, only: BndType
+    use RivModule, only: RivType
     class(BndType), pointer :: packobj
     ! -- dummy
     real(DP), intent(in), dimension(:) :: hnew
@@ -674,39 +675,42 @@ subroutine buy_cf_riv(packobj, hnew, dense, elev, denseref, locelev, &
 ! ------------------------------------------------------------------------------
     !
     ! -- Process density terms for each RIV
-    do n = 1, packobj%nbound
-      node = packobj%nodelist(n)
-      if (packobj%ibound(node) <= 0) cycle
-      !
-      ! -- density
-      denseriv = get_bnd_density(n, locdense, locconc, denseref, &
-                                 drhodc, crhoref, ctemp, packobj%auxvar)
-      !
-      ! -- elevation
-      elevriv = elev(node)
-      if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
-      !
-      ! -- boundary head and conductance
-      hriv = packobj%bound(1, n)
-      cond = packobj%bound(2, n)
-      rbot = packobj%bound(3, n)
-      !
-      ! -- calculate and add terms depending on whether head is above rbot
-      if (hnew(node) > rbot) then
+    select type (packobj)
+    type is (RivType)
+      do n = 1, packobj%nbound
+        node = packobj%nodelist(n)
+        if (packobj%ibound(node) <= 0) cycle
         !
-        ! --calculate HCOF and RHS terms, similar to GHB in this case
-        call calc_ghb_hcof_rhs_terms(denseref, denseriv, dense(node), &
-                                     elevriv, elev(node), hriv, hnew(node), &
-                                     cond, iform, rhsterm, hcofterm)
-      else
-        hcofterm = DZERO
-        rhsterm = cond * (denseriv / denseref - DONE) * (hriv - rbot)
-      end if
-      !
-      ! -- Add terms to package hcof and rhs accumulators
-      packobj%hcof(n) = packobj%hcof(n) + hcofterm
-      packobj%rhs(n) = packobj%rhs(n) - rhsterm
-    end do
+        ! -- density
+        denseriv = get_bnd_density(n, locdense, locconc, denseref, &
+                                   drhodc, crhoref, ctemp, packobj%auxvar)
+        !
+        ! -- elevation
+        elevriv = elev(node)
+        if (locelev > 0) elevriv = packobj%auxvar(locelev, n)
+        !
+        ! -- boundary head and conductance
+        hriv = packobj%stage(n)
+        cond = packobj%cond(n)
+        rbot = packobj%rbot(n)
+        !
+        ! -- calculate and add terms depending on whether head is above rbot
+        if (hnew(node) > rbot) then
+          !
+          ! --calculate HCOF and RHS terms, similar to GHB in this case
+          call calc_ghb_hcof_rhs_terms(denseref, denseriv, dense(node), &
+                                       elevriv, elev(node), hriv, hnew(node), &
+                                       cond, iform, rhsterm, hcofterm)
+        else
+          hcofterm = DZERO
+          rhsterm = cond * (denseriv / denseref - DONE) * (hriv - rbot)
+        end if
+        !
+        ! -- Add terms to package hcof and rhs accumulators
+        packobj%hcof(n) = packobj%hcof(n) + hcofterm
+        packobj%rhs(n) = packobj%rhs(n) - rhsterm
+      end do
+    end select
     !
     ! -- Return
     return