feat(mst): add capability to save sorbate concentrations to binary file

autotest/test_gwt_mst02.py

@@ -1,5 +1,8 @@
 """
-Test the GWT Sorption (RCT) Package by running a ...
+Test the GWT Sorption (RCT) Package by running a simple 2-cell test with
+mass injected into the first cell at a rate of 1.0 unit/time.  Compare the
+simulated aqueous and sorbate concentrations with a known solution.  Problem
+uses 10 time steps to discretize a 1.0 time unit period.
 """
 
 import os
@@ -189,12 +192,13 @@ def build_models(idx, test):
         porosity=sy,
         sorption="linear",
         bulk_density=1.0,
+        sorbate_filerecord=f"{gwtname}.mst.csrb",
         distcoef=distcoef[idx],
     )
 
     # mass loading source
     srcdict = {0: [[(0, 0, 0), 1.0]]}
-    cnc = flopy.mf6.ModflowGwtsrc(
+    src = flopy.mf6.ModflowGwtsrc(
         gwt, stress_period_data=srcdict, save_flows=False, pname="SRC-1"
     )
 
@@ -229,6 +233,7 @@ def check_output(idx, test):
     name = cases[idx]
     gwtname = "gwt_" + name
 
+    # Check aqueous concentrations
     fpth = os.path.join(test.workspace, f"{gwtname}.ucn")
     try:
         cobj = flopy.utils.HeadFile(
@@ -237,10 +242,23 @@ def check_output(idx, test):
         ts = cobj.get_ts([(0, 0, 0), (0, 0, 1)])
     except:
         assert False, f'could not load data from "{fpth}"'
+    assert np.allclose(
+        ts, tsanswers[idx]
+    ), "simulated concentrations do not match with known solution."
 
-    # Check concentrations
-    assert np.allclose(ts, tsanswers[idx]), (
-        "simulated concentrations do not " "match with known solution."
+    # Check sorbate concentrations
+    fpth = os.path.join(test.workspace, f"{gwtname}.mst.csrb")
+    try:
+        cobj = flopy.utils.HeadFile(fpth, precision="double", text="SORBATE")
+        ts = cobj.get_ts([(0, 0, 0), (0, 0, 1)])
+    except:
+        assert False, f'could not load data from "{fpth}"'
+    d = distcoef[idx]
+    tsa_csrb = tsanswers[idx]
+    tsa_csrb[:, 1:] *= d
+    assert np.allclose(ts, tsa_csrb), (
+        "Sorbate concentrations do not match with known solution.",
+        ts,
     )
 
     # Check budget file

doc/ReleaseNotes/develop.tex

@@ -6,7 +6,7 @@
 	\underline{NEW FUNCTIONALITY}
 	\begin{itemize}
 		\item A new adaptive time stepping (ATS) capability was added to the Advection (ADV) Package of the Groundwater Transport (GWT) Model.  A new input option, called ATS\_PERCEL, specifies the fractional cell distance that a particle of water can travel within one time step.  When ATS\_PERCEL is specified by the user, and the ATS utility is activated in the TDIS Package, the ADV Package will calculate the largest time step that will meet this fractional cell distance constraint, and will submit this time step to the ATS utility.  This option may improve time stepping for solute transport models and for variable-density flow and transport models by allowing step lengths to be calculated as a function of the flow system rather than being specified as input by the user.
-	%	\item xxx
+		\item Added the capability to write sorbate concentrations to a binary output file.  A new SORBATE option is now available in the Mass Storage and Transfer (MST) Package of the GWT Model to provide the name of the binary output file for sorbate concentrations. Sorbate concentrations will be written to the binary output file whenever concentrations for the GWT model are saved, as determined by the GWT Output Control option.
 	%	\item xxx
 	\end{itemize}
 

doc/mf6io/mf6ivar/dfn/gwt-mst.dfn

@@ -33,6 +33,50 @@ optional true
 longname activate sorption
 description is a text keyword to indicate that sorption will be activated.  Valid sorption options include LINEAR, FREUNDLICH, and LANGMUIR.  Use of this keyword requires that BULK\_DENSITY and DISTCOEF are specified in the GRIDDATA block.  If sorption is specified as FREUNDLICH or LANGMUIR then SP2 is also required in the GRIDDATA block.
 
+block options
+name sorbate_filerecord
+type record sorbate fileout sorbatefile
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
+block options
+name sorbate
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname sorbate keyword
+description keyword to specify that record corresponds to sorbate concentration.
+
+block options
+name fileout
+type keyword
+shape
+in_record true
+reader urword
+tagged true
+optional false
+longname file keyword
+description keyword to specify that an output filename is expected next.
+
+block options
+name sorbatefile
+type string
+preserve_case true
+shape
+in_record true
+reader urword
+tagged false
+optional false
+longname file keyword
+description name of the output file to write sorbate concentration information.  Sorbate concentrations will be written whenever aqueous concentrations are saved, as determined by settings in the Output Control option.
+
 # --------------------- gwt mst griddata ---------------------
 
 block griddata

src/Model/GroundWaterEnergy/gwe.f90

@@ -50,7 +50,7 @@ module GweModule
     procedure :: model_cc => gwe_cc
     procedure :: model_cq => gwe_cq
     procedure :: model_bd => gwe_bd
-    procedure :: model_ot => gwe_ot
+    procedure :: tsp_ot_flow => gwe_ot_flow
     procedure :: model_da => gwe_da
     procedure :: model_bdentry => gwe_bdentry
     procedure :: allocate_scalars
@@ -568,32 +568,22 @@ subroutine gwe_bd(this, icnvg, isuppress_output)
     return
   end subroutine gwe_bd
 
-  !> @brief GWE Model Output
+  !> @brief GWE model output routine
   !!
-  !! This subroutine calls the parent class output routine.
+  !! Save and print flows
   !<
-  subroutine gwe_ot(this)
-    ! -- dummy
+  subroutine gwe_ot_flow(this, icbcfl, ibudfl, icbcun)
+    ! dummy
     class(GweModelType) :: this
-    ! -- local
-    integer(I4B) :: icbcfl
-    integer(I4B) :: icbcun
-    ! -- formats
-    !
-    ! -- Initialize
-    icbcfl = 0
-    !
-    ! -- Because est belongs to gwe, call est_ot_flow directly (and not from parent)
-    if (this%oc%oc_save('BUDGET')) icbcfl = 1
-    icbcun = this%oc%oc_save_unit('BUDGET')
+    integer(I4B), intent(in) :: icbcfl
+    integer(I4B), intent(in) :: ibudfl
+    integer(I4B), intent(in) :: icbcun
+    ! local
+
     if (this%inest > 0) call this%est%est_ot_flow(icbcfl, icbcun)
-    !
-    ! -- Call parent class _ot routines.
-    call this%tsp_ot(this%inest)
-    !
-    ! -- Return
-    return
-  end subroutine gwe_ot
+    call this%TransportModelType%tsp_ot_flow(icbcfl, ibudfl, icbcun)
+
+  end subroutine gwe_ot_flow
 
   !> @brief Deallocate
   !!

src/Model/GroundWaterTransport/gwt-mst.f90

@@ -10,7 +10,8 @@
 module GwtMstModule
 
   use KindModule, only: DP, I4B
-  use ConstantsModule, only: DONE, DZERO, DTWO, DHALF, LENBUDTXT
+  use ConstantsModule, only: DONE, DZERO, DTWO, DHALF, LENBUDTXT, MAXCHARLEN, &
+                             MNORMAL, LINELENGTH, DHNOFLO
   use SimVariablesModule, only: errmsg, warnmsg
   use SimModule, only: store_error, count_errors, &
                        store_warning
@@ -52,11 +53,13 @@ module GwtMstModule
     !
     ! -- sorption
     integer(I4B), pointer :: isrb => null() !< sorption active flag (0:off, 1:linear, 2:freundlich, 3:langmuir)
+    integer(I4B), pointer :: ioutsorbate => null() !< unit number for sorbate concentration output
     real(DP), dimension(:), pointer, contiguous :: bulk_density => null() !< bulk density of mobile domain; mass of mobile domain solid per aquifer volume
     real(DP), dimension(:), pointer, contiguous :: distcoef => null() !< kd distribution coefficient
     real(DP), dimension(:), pointer, contiguous :: sp2 => null() !< second sorption parameter
     real(DP), dimension(:), pointer, contiguous :: ratesrb => null() !< rate of sorption
     real(DP), dimension(:), pointer, contiguous :: ratedcys => null() !< rate of sorbed mass decay
+    real(DP), dimension(:), pointer, contiguous :: csrb => null() !< sorbate concentration
     !
     ! -- misc
     integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound
@@ -75,8 +78,10 @@ module GwtMstModule
     procedure :: mst_cq_dcy
     procedure :: mst_cq_srb
     procedure :: mst_cq_dcy_srb
+    procedure :: mst_calc_csrb
     procedure :: mst_bd
     procedure :: mst_ot_flow
+    procedure :: mst_ot_dv
     procedure :: mst_da
     procedure :: allocate_scalars
     procedure :: addto_volfracim
@@ -590,6 +595,11 @@ subroutine mst_cq(this, nodes, cnew, cold, flowja)
       call this%mst_cq_dcy_srb(nodes, cnew, cold, flowja)
     end if
     !
+    ! -- calculate csrb
+    if (this%isrb /= 0) then
+      call this%mst_calc_csrb(cnew)
+    end if
+    !
     ! -- Return
     return
   end subroutine mst_cq
@@ -867,10 +877,36 @@ subroutine mst_cq_dcy_srb(this, nodes, cnew, cold, flowja)
     return
   end subroutine mst_cq_dcy_srb
 
+  !> @ brief Calculate sorbed concentration
+  !<
+  subroutine mst_calc_csrb(this, cnew)
+    ! -- dummy
+    class(GwtMstType) :: this !< GwtMstType object
+    real(DP), intent(in), dimension(:) :: cnew !< concentration at end of this time step
+    ! -- local
+    integer(I4B) :: n
+    real(DP) :: distcoef
+    real(DP) :: csrb
+
+    ! Calculate sorbed concentration
+    do n = 1, size(cnew)
+      csrb = DZERO
+      if (this%ibound(n) > 0) then
+        distcoef = this%distcoef(n)
+        if (this%isrb == 1) then
+          csrb = cnew(n) * distcoef
+        else if (this%isrb == 2) then
+          csrb = get_freundlich_conc(cnew(n), distcoef, this%sp2(n))
+        else if (this%isrb == 3) then
+          csrb = get_langmuir_conc(cnew(n), distcoef, this%sp2(n))
+        end if
+      end if
+      this%csrb(n) = csrb
+    end do
+
+  end subroutine mst_calc_csrb
+
   !> @ brief Calculate budget terms for package
-  !!
-  !!  Method to calculate budget terms for the package.
-  !!
   !<
   subroutine mst_bd(this, isuppress_output, model_budget)
     ! -- modules
@@ -974,6 +1010,42 @@ subroutine mst_ot_flow(this, icbcfl, icbcun)
     return
   end subroutine mst_ot_flow
 
+  !> @brief Save sorbate concentration array to binary file
+  !<
+  subroutine mst_ot_dv(this, idvsave)
+    ! -- dummy
+    class(GwtMstType) :: this
+    integer(I4B), intent(in) :: idvsave
+    ! -- local
+    character(len=1) :: cdatafmp = ' ', editdesc = ' '
+    integer(I4B) :: ibinun
+    integer(I4B) :: iprint
+    integer(I4B) :: nvaluesp
+    integer(I4B) :: nwidthp
+    real(DP) :: dinact
+
+    ! Set unit number for sorbate output
+    if (this%ioutsorbate /= 0) then
+      ibinun = 1
+    else
+      ibinun = 0
+    end if
+    if (idvsave == 0) ibinun = 0
+
+    ! save sorbate concentration array
+    if (ibinun /= 0) then
+      iprint = 0
+      dinact = DHNOFLO
+      if (this%ioutsorbate /= 0) then
+        ibinun = this%ioutsorbate
+        call this%dis%record_array(this%csrb, this%iout, iprint, ibinun, &
+                                   '         SORBATE', cdatafmp, nvaluesp, &
+                                   nwidthp, editdesc, dinact)
+      end if
+    end if
+
+  end subroutine mst_ot_dv
+
   !> @ brief Deallocate
   !!
   !!  Method to deallocate memory for the package.
@@ -998,10 +1070,12 @@ subroutine mst_da(this)
       call mem_deallocate(this%decaylast)
       call mem_deallocate(this%decayslast)
       call mem_deallocate(this%isrb)
+      call mem_deallocate(this%ioutsorbate)
       call mem_deallocate(this%bulk_density)
       call mem_deallocate(this%distcoef)
       call mem_deallocate(this%sp2)
       call mem_deallocate(this%ratesrb)
+      call mem_deallocate(this%csrb)
       call mem_deallocate(this%ratedcys)
       this%ibound => null()
       this%fmi => null()
@@ -1033,10 +1107,12 @@ subroutine allocate_scalars(this)
     !
     ! -- Allocate
     call mem_allocate(this%isrb, 'ISRB', this%memoryPath)
+    call mem_allocate(this%ioutsorbate, 'IOUTSORBATE', this%memoryPath)
     call mem_allocate(this%idcy, 'IDCY', this%memoryPath)
     !
     ! -- Initialize
     this%isrb = 0
+    this%ioutsorbate = 0
     this%idcy = 0
     !
     ! -- Return
@@ -1093,10 +1169,12 @@ subroutine allocate_arrays(this, nodes)
       call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
       call mem_allocate(this%distcoef, 1, 'DISTCOEF', this%memoryPath)
       call mem_allocate(this%ratesrb, 1, 'RATESRB', this%memoryPath)
+      call mem_allocate(this%csrb, 1, 'CSRB', this%memoryPath)
     else
       call mem_allocate(this%bulk_density, nodes, 'BULK_DENSITY', this%memoryPath)
       call mem_allocate(this%distcoef, nodes, 'DISTCOEF', this%memoryPath)
       call mem_allocate(this%ratesrb, nodes, 'RATESRB', this%memoryPath)
+      call mem_allocate(this%csrb, nodes, 'CSRB', this%memoryPath)
       if (this%isrb == 1) then
         call mem_allocate(this%sp2, 1, 'SP2', this%memoryPath)
       else
@@ -1120,6 +1198,7 @@ subroutine allocate_arrays(this, nodes)
       this%bulk_density(n) = DZERO
       this%distcoef(n) = DZERO
       this%ratesrb(n) = DZERO
+      this%csrb(n) = DZERO
     end do
     do n = 1, size(this%sp2)
       this%sp2(n) = DZERO
@@ -1140,11 +1219,13 @@ end subroutine allocate_arrays
   !<
   subroutine read_options(this)
     ! -- modules
-    use ConstantsModule, only: LINELENGTH
+    use OpenSpecModule, only: access, form
+    use InputOutputModule, only: getunit, openfile
     ! -- dummy
     class(GwtMstType) :: this !< GwtMstType object
     ! -- local
     character(len=LINELENGTH) :: keyword, keyword2
+    character(len=MAXCHARLEN) :: fname
     integer(I4B) :: ierr
     logical :: isfound, endOfBlock
     ! -- formats
@@ -1161,6 +1242,9 @@ subroutine read_options(this)
       &"(4x,'FIRST-ORDER DECAY IS ACTIVE. ')"
     character(len=*), parameter :: fmtidcy2 = &
       &"(4x,'ZERO-ORDER DECAY IS ACTIVE. ')"
+    character(len=*), parameter :: fmtfileout = &
+      "(4x,'MST ',1x,a,1x,' WILL BE SAVED TO FILE: ',a,/4x,&
+      &'OPENED ON UNIT: ',I7)"
     !
     ! -- get options block
     call this%parser%GetBlock('OPTIONS', isfound, ierr, blockRequired=.false., &
@@ -1197,6 +1281,20 @@ subroutine read_options(this)
         case ('ZERO_ORDER_DECAY')
           this%idcy = 2
           write (this%iout, fmtidcy2)
+        case ('SORBATE')
+          call this%parser%GetStringCaps(keyword)
+          if (keyword == 'FILEOUT') then
+            call this%parser%GetString(fname)
+            this%ioutsorbate = getunit()
+            call openfile(this%ioutsorbate, this%iout, fname, 'DATA(BINARY)', &
+                          form, access, 'REPLACE', mode_opt=MNORMAL)
+            write (this%iout, fmtfileout) &
+              'SORBATE', fname, this%ioutsorbate
+          else
+            errmsg = 'Optional SORBATE keyword must be '// &
+                     'followed by FILEOUT.'
+            call store_error(errmsg)
+          end if
         case default
           write (errmsg, '(a,a)') 'Unknown MST option: ', trim(keyword)
           call store_error(errmsg)