diff --git a/autotest/test_gwt_mst02.py b/autotest/test_gwt_mst02.py index 3dc96596ce0..c4898b871f3 100644 --- a/autotest/test_gwt_mst02.py +++ b/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 diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index e9cd8590a3f..d6a0e3764a1 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/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} diff --git a/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn b/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn index 06e5169f816..dbaaf5c6fe9 100644 --- a/doc/mf6io/mf6ivar/dfn/gwt-mst.dfn +++ b/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 diff --git a/src/Model/GroundWaterEnergy/gwe.f90 b/src/Model/GroundWaterEnergy/gwe.f90 index 947584b24c8..92b4efc530d 100644 --- a/src/Model/GroundWaterEnergy/gwe.f90 +++ b/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 !! diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index bc58193778a..d20fb7d15ed 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/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) diff --git a/src/Model/GroundWaterTransport/gwt.f90 b/src/Model/GroundWaterTransport/gwt.f90 index 9a49e077276..9c8f012b017 100644 --- a/src/Model/GroundWaterTransport/gwt.f90 +++ b/src/Model/GroundWaterTransport/gwt.f90 @@ -54,7 +54,8 @@ module GwtModule procedure :: model_cc => gwt_cc procedure :: model_cq => gwt_cq procedure :: model_bd => gwt_bd - procedure :: model_ot => gwt_ot + procedure :: tsp_ot_flow => gwt_ot_flow + procedure :: tsp_ot_dv => gwt_ot_dv procedure :: model_da => gwt_da procedure :: model_bdentry => gwt_bdentry procedure :: allocate_scalars @@ -603,32 +604,41 @@ subroutine gwt_bd(this, icnvg, isuppress_output) return end subroutine gwt_bd - !> @brief Print and/or save model output + !> @brief GWT model output routine !! - !! Call the parent class output routine + !! Save and print flows !< - subroutine gwt_ot(this) - ! -- dummy + subroutine gwt_ot_flow(this, icbcfl, ibudfl, icbcun) + ! dummy class(GwtModelType) :: this - ! -- local - integer(I4B) :: icbcfl - integer(I4B) :: icbcun - ! - ! - ! -- Initialize - icbcfl = 0 - ! - ! -- Because mst belongs to gwt, call mst_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 + + ! call GWT flow output routines if (this%inmst > 0) call this%mst%mst_ot_flow(icbcfl, icbcun) - ! - ! -- Call parent class _ot routines. - call this%tsp_ot(this%inmst) - ! - ! -- Return - return - end subroutine gwt_ot + + ! call general transport model flow output routines + call this%TransportModelType%tsp_ot_flow(icbcfl, ibudfl, icbcun) + + end subroutine gwt_ot_flow + + !> @brief GWT model dependent variable output + !< + subroutine gwt_ot_dv(this, idvsave, idvprint, ipflag) + class(GwtModelType) :: this + integer(I4B), intent(in) :: idvsave + integer(I4B), intent(in) :: idvprint + integer(I4B), intent(inout) :: ipflag + + ! call GWT dependent variable output routines + if (this%inmst > 0) call this%mst%mst_ot_dv(idvsave) + + ! call general transport model dependent variable output routines + call this%TransportModelType%tsp_ot_dv(idvsave, idvprint, ipflag) + + end subroutine gwt_ot_dv !> @brief Deallocate !! diff --git a/src/Model/TransportModel/tsp.f90 b/src/Model/TransportModel/tsp.f90 index 7020c8b8d78..c0a172f82cf 100644 --- a/src/Model/TransportModel/tsp.f90 +++ b/src/Model/TransportModel/tsp.f90 @@ -68,15 +68,15 @@ module TransportModelModule procedure, public :: tsp_cc procedure, public :: tsp_cq procedure, public :: tsp_bd - procedure, public :: tsp_ot + procedure, public :: model_ot => tsp_ot + procedure, public :: tsp_ot_flow + procedure, public :: tsp_ot_dv procedure, public :: allocate_tsp_scalars procedure, public :: set_tsp_labels procedure, public :: ftype_check ! -- private procedure, private :: tsp_ot_obs - procedure, private :: tsp_ot_flow procedure, private :: tsp_ot_flowja - procedure, private :: tsp_ot_dv procedure, private :: tsp_ot_bdsummary procedure, private :: create_tsp_packages procedure, private :: log_namfile_options @@ -304,12 +304,11 @@ end subroutine tsp_bd !! !! Generalized transport model output !< - subroutine tsp_ot(this, inmst) + subroutine tsp_ot(this) ! -- modules use TdisModule, only: kstp, kper, tdis_ot, endofperiod ! -- dummy class(TransportModelType) :: this - integer(I4B), intent(in) :: inmst ! -- local integer(I4B) :: idvsave integer(I4B) :: idvprint @@ -343,7 +342,7 @@ subroutine tsp_ot(this, inmst) call this%tsp_ot_obs() ! ! -- Save and print flows - call this%tsp_ot_flow(icbcfl, ibudfl, icbcun, inmst) + call this%tsp_ot_flow(icbcfl, ibudfl, icbcun) ! ! -- Save and print dependent variables call this%tsp_ot_dv(idvsave, idvprint, ipflag) @@ -389,13 +388,12 @@ end subroutine tsp_ot_obs !! !! Save and print flows !< - subroutine tsp_ot_flow(this, icbcfl, ibudfl, icbcun, inmst) + subroutine tsp_ot_flow(this, icbcfl, ibudfl, icbcun) ! -- dummy class(TransportModelType) :: this integer(I4B), intent(in) :: icbcfl integer(I4B), intent(in) :: ibudfl integer(I4B), intent(in) :: icbcun - integer(I4B), intent(in) :: inmst ! -- local class(BndType), pointer :: packobj integer(I4B) :: ip @@ -479,7 +477,7 @@ subroutine tsp_ot_flowja(this, nja, flowja, icbcfl, icbcun) return end subroutine tsp_ot_flowja - !> @brief Generalized tranpsort model output routine + !> @brief Generalized transport model output routine !! !! Loop through attached packages saving and printing dependent variables !<