feat(idm): update RCH package for IDM (MODFLOW-USGS#1380)

* update RCH package for IDM

* rebuild makefiles

* some rch package cleanup

* add write_list call, additional cleanup

---------

Co-authored-by: mjreno <mreno@IGSAAA071L00066.gs.doi.net>

autotest/test_gwf_libmf6_rch01.py

@@ -177,7 +177,7 @@ def api_func(exe, idx, model_ws=None):
     max_iter = mf6.get_value(mxit_tag)
 
     # get copy of recharge array
-    rch_tag = mf6.get_var_address("BOUND", name, rch_pname)
+    rch_tag = mf6.get_var_address("RECHARGE", name, rch_pname)
     new_recharge = mf6.get_value(rch_tag)
 
     # model time loop
@@ -193,7 +193,7 @@ def api_func(exe, idx, model_ws=None):
         mf6.prepare_solve()
 
         # update recharge
-        new_recharge[:, 0] = rch_spd[idx] * area
+        new_recharge[:] = rch_spd[idx]
         mf6.set_value(rch_tag, new_recharge)
 
         while kiter < max_iter:

autotest/test_gwf_libmf6_rch02.py

@@ -172,7 +172,7 @@ def run_perturbation(mf6, max_iter, recharge, tag, rch):
     kiter = 0
     while kiter < max_iter:
         # update recharge
-        recharge[:, 0] = rch * area
+        recharge[:] = rch
         mf6.set_value(tag, recharge)
         # solve with updated well rate
         has_converged = mf6.solve()
@@ -222,7 +222,7 @@ def api_func(exe, idx, model_ws=None):
     max_iter = mf6.get_value(mxit_tag)
 
     # get copy of recharge array
-    rch_tag = mf6.get_var_address("BOUND", name, rch_pname)
+    rch_tag = mf6.get_var_address("RECHARGE", name, rch_pname)
     new_recharge = mf6.get_value(rch_tag).copy()
 
     # determine initial recharge value

doc/mf6io/mf6ivar/dfn/gwf-rch.dfn

@@ -1,5 +1,6 @@
 # --------------------- gwf rch options ---------------------
 # flopy multi-package
+# modflow6 aux-sfac-param recharge
 
 block options
 name fixed_cell
@@ -44,6 +45,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'recharge'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -52,6 +54,7 @@ reader urword
 optional true
 longname print recharge rates to listing file
 description REPLACE print_flows {'{#1}': 'recharge'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -60,6 +63,7 @@ reader urword
 optional true
 longname save recharge to budget file
 description REPLACE save_flows {'{#1}': 'recharge'}
+mf6internal ipakcb
 
 block options
 name ts_filerecord
@@ -169,6 +173,7 @@ shape (maxbound)
 reader urword
 longname
 description
+mf6internal spd
 
 block period
 name cellid
@@ -202,6 +207,7 @@ optional true
 time_series true
 longname auxiliary variables
 description REPLACE aux {'{#1}': 'recharge'}
+mf6internal auxvar
 
 block period
 name boundname

doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn

@@ -1,5 +1,6 @@
 # --------------------- gwf rcha options ---------------------
 # flopy multi-package
+# modflow6 aux-sfac-param recharge
 
 block options
 name readasarrays
@@ -45,6 +46,7 @@ reader urword
 optional true
 longname print input to listing file
 description REPLACE print_input {'{#1}': 'recharge'}
+mf6internal iprpak
 
 block options
 name print_flows
@@ -53,6 +55,7 @@ reader urword
 optional true
 longname print recharge rates to listing file
 description REPLACE print_flows {'{#1}': 'recharge'}
+mf6internal iprflow
 
 block options
 name save_flows
@@ -61,6 +64,7 @@ reader urword
 optional true
 longname save CHD flows to budget file
 description REPLACE save_flows {'{#1}': 'recharge'}
+mf6internal ipakcb
 
 block options
 name tas_filerecord
@@ -168,6 +172,7 @@ name recharge
 type double precision
 shape (ncol*nrow; ncpl)
 reader readarray
+time_series true
 longname recharge rate
 description is the recharge flux rate ($LT^{-1}$).  This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the "Using Time-Array Series in a Package" section).
 default_value 1.e-3
@@ -177,6 +182,8 @@ name aux
 type double precision
 shape (ncol*nrow; ncpl)
 reader readarray
+time_series true
 optional true
 longname auxiliary variable iaux
 description is an array of values for auxiliary variable aux(iaux), where iaux is a value from 1 to naux, and aux(iaux) must be listed as part of the auxiliary variables.  A separate array can be specified for each auxiliary variable.  If an array is not specified for an auxiliary variable, then a value of zero is assigned.  If the value specified here for the auxiliary variable is the same as auxmultname, then the recharge array will be multiplied by this array.
+mf6internal auxvar

make/makefile

@@ -104,6 +104,8 @@ $(OBJDIR)/gwt1disu1idm.o \
 $(OBJDIR)/gwt1dis1idm.o \
 $(OBJDIR)/gwf3wel8idm.o \
 $(OBJDIR)/gwf3riv8idm.o \
+$(OBJDIR)/gwf3rch8idm.o \
+$(OBJDIR)/gwf3rcha8idm.o \
 $(OBJDIR)/gwf3npf8idm.o \
 $(OBJDIR)/gwf3idm.o \
 $(OBJDIR)/gwf3ghb8idm.o \

msvs/mf6core.vfproj

@@ -146,6 +146,8 @@
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3obs8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3oc8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rch8idm.f90"/>
+		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3rcha8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3riv8idm.f90"/>
 		<File RelativePath="..\src\Model\GroundWaterFlow\gwf3sfr8.f90"/>

src/Model/GroundWaterFlow/gwf3.f90

@@ -1305,7 +1305,8 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
       call ghb_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
                       pakname, mempath)
     case ('RCH6')
-      call rch_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
+      call rch_create(packobj, ipakid, ipaknum, inunit, iout, this%name, &
+                      pakname, mempath)
     case ('EVT6')
       call evt_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname)
     case ('MAW6')

src/Model/GroundWaterFlow/gwf3dis8.f90

@@ -36,6 +36,7 @@ module GwfDisModule
     procedure, public :: read_layer_array
     procedure, public :: record_srcdst_list_header
     procedure, public :: nlarray_to_nodelist
+    procedure, public :: nlarray_to_nodelist2
     ! -- helper functions
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx3
@@ -1747,4 +1748,91 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, &
     ! -- return
   end subroutine nlarray_to_nodelist
 
+  subroutine nlarray_to_nodelist2(this, darray, nodelist, maxbnd, nbound, aname)
+! ******************************************************************************
+! nlarray_to_nodelist -- Read an integer array into nodelist. For structured
+!                        model, integer array is layer number; for unstructured
+!                        model, integer array is node number.
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use InputOutputModule, only: get_node
+    use ConstantsModule, only: LINELENGTH
+    ! -- dummy
+    class(GwfDisType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    ! -- local
+    integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
+! ------------------------------------------------------------------------------
+    !
+    ! -- set variables
+    nlay = this%mshape(1)
+    nrow = this%mshape(2)
+    ncol = this%mshape(3)
+    !
+    if (this%ndim > 1) then
+      !
+      nval = ncol * nrow
+      !
+      ! -- Copy array into nodelist
+      ipos = 1
+      ierr = 0
+      do ir = 1, nrow
+        do ic = 1, ncol
+          nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
+          il = darray(nodeu)
+          if (il < 1 .or. il > nlay) then
+            write (errmsg, '(a,1x,i0)') 'Invalid layer number:', il
+            call store_error(errmsg, terminate=.TRUE.)
+          end if
+          nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
+          noder = this%get_nodenumber(nodeu, 0)
+          if (ipos > maxbnd) then
+            ierr = ipos
+          else
+            nodelist(ipos) = noder
+          end if
+          ipos = ipos + 1
+        end do
+      end do
+      !
+      ! -- Check for errors
+      nbound = ipos - 1
+      if (ierr > 0) then
+        write (errmsg, '(a,1x,i0)') &
+          'MAXBOUND dimension is too small.'// &
+          'INCREASE MAXBOUND TO:', ierr
+        call store_error(errmsg, terminate=.TRUE.)
+      end if
+      !
+      ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
+      if (nbound < maxbnd) then
+        do ipos = nbound + 1, maxbnd
+          nodelist(ipos) = 0
+        end do
+      end if
+      !
+    else
+      !
+      ! -- For unstructured, read nodelist directly, then check node numbers
+      nodelist = darray
+      do noder = 1, maxbnd
+        if (noder < 1 .or. noder > this%nodes) then
+          write (errmsg, '(a,1x,i0)') 'Invalid node number:', noder
+          call store_error(errmsg, terminate=.TRUE.)
+        end if
+      end do
+      nbound = maxbnd
+      !
+    end if
+    !
+    ! -- return
+  end subroutine nlarray_to_nodelist2
+
 end module GwfDisModule

src/Model/GroundWaterFlow/gwf3disv8.f90

@@ -38,6 +38,7 @@ module GwfDisvModule
     procedure, public :: read_layer_array
     procedure, public :: record_srcdst_list_header
     procedure, public :: nlarray_to_nodelist
+    procedure, public :: nlarray_to_nodelist2
     ! -- helper functions
     procedure :: get_nodenumber_idx1
     procedure :: get_nodenumber_idx2
@@ -1981,4 +1982,74 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, &
     ! -- return
   end subroutine nlarray_to_nodelist
 
+  subroutine nlarray_to_nodelist2(this, darray, nodelist, maxbnd, nbound, aname)
+! ******************************************************************************
+! nlarray_to_nodelist -- Read an integer array into nodelist. For structured
+!                        model, integer array is layer number; for unstructured
+!                        model, integer array is node number.
+! ******************************************************************************
+!
+!    SPECIFICATIONS:
+! ------------------------------------------------------------------------------
+    ! -- modules
+    use InputOutputModule, only: get_node
+    ! -- dummy
+    class(GwfDisvType) :: this
+    integer(I4B), intent(in) :: maxbnd
+    integer(I4B), dimension(:), pointer, contiguous :: darray
+    integer(I4B), dimension(maxbnd), intent(inout) :: nodelist
+    integer(I4B), intent(inout) :: nbound
+    character(len=*), intent(in) :: aname
+    ! -- local
+    integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr
+! ------------------------------------------------------------------------------
+    !
+    ! -- set variables
+    nlay = this%mshape(1)
+    nrow = 1
+    ncol = this%mshape(2)
+    !
+    nval = ncol * nrow
+    !
+    ! -- Copy array into nodelist
+    ipos = 1
+    ierr = 0
+    do ir = 1, nrow
+      do ic = 1, ncol
+        nodeu = get_node(1, ir, ic, nlay, nrow, ncol)
+        il = darray(nodeu)
+        if (il < 1 .or. il > nlay) then
+          write (errmsg, '(a,i0,a)') &
+            'Invalid layer number (', il, ').'
+          call store_error(errmsg, terminate=.TRUE.)
+        end if
+        nodeu = get_node(il, ir, ic, nlay, nrow, ncol)
+        noder = this%get_nodenumber(nodeu, 0)
+        if (ipos > maxbnd) then
+          ierr = ipos
+        else
+          nodelist(ipos) = noder
+        end if
+        ipos = ipos + 1
+      end do
+    end do
+    !
+    ! -- Check for errors
+    nbound = ipos - 1
+    if (ierr > 0) then
+      write (errmsg, '(a,i0,a)') &
+        'MAXBOUND dimension is too small. Increase MAXBOUND to ', ierr, '.'
+      call store_error(errmsg, terminate=.TRUE.)
+    end if
+    !
+    ! -- If nbound < maxbnd, then initialize nodelist to zero in this range
+    if (nbound < maxbnd) then
+      do ipos = nbound + 1, maxbnd
+        nodelist(ipos) = 0
+      end do
+    end if
+    !
+    ! -- return
+  end subroutine nlarray_to_nodelist2
+
 end module GwfDisvModule