From 20ae3f12b99cd77e2b78a6fc8e008f7c88861cb6 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 20 Jan 2025 19:00:02 -0800 Subject: [PATCH 01/18] fix(gwe-est): add support for energy decay in the solid phase --- autotest/test_gwe_decay01.py | 177 +++++++++++ doc/ReleaseNotes/develop.tex | 1 + doc/mf6io/mf6ivar/dfn/gwe-est.dfn | 28 +- src/Model/GroundWaterEnergy/gwe-est.f90 | 323 ++++++++++++++++----- src/Model/GroundWaterTransport/gwt-mst.f90 | 1 - 5 files changed, 457 insertions(+), 73 deletions(-) create mode 100644 autotest/test_gwe_decay01.py diff --git a/autotest/test_gwe_decay01.py b/autotest/test_gwe_decay01.py new file mode 100644 index 00000000000..7999b6b1e02 --- /dev/null +++ b/autotest/test_gwe_decay01.py @@ -0,0 +1,177 @@ +""" +Test problem for decay of energy in EST package of GWE. Uses a single-cell +model. Test contributed by Cas Neyens. +""" + +# Imports + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +# Base simulation and model name and workspace + +cases = ["decay-aqe", "decay-sld", "decay-both"] + +# Model units +length_units = "meters" +time_units = "seconds" + +rho_w = 1000 # kg/m^3 +rho_s = 2500 # kg/m^3 +n = 0.2 # - +gamma_w = -1000 # J/s/m^3, arbitrary value for zero-order aqueous heat production +gamma_s = -0.1 # J/s/kg, arbitrary value for zero-order solid heat production +c_w = 4000 # J/kg/degC +c_s = 1000 # J/kg/decC +T0 = 0 # degC + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1 # m +delc = 1 # m +top = 1 # m +botm = 0 # m + +perlen = 86400 # s +nstp = 20 + +parameters = { + # aqueous + "decay-aqe": { + "zero_order_decay_water": True, + "zero_order_decay_solid": False, + "decay_water": gamma_w, + "decay_solid": gamma_s, + }, + # solid + "decay-sld": { + "zero_order_decay_water": False, + "zero_order_decay_solid": True, + "decay_water": gamma_w, + "decay_solid": gamma_s, + }, + # combined + "decay-both": { + "zero_order_decay_water": True, + "zero_order_decay_solid": True, + "decay_water": gamma_w, + "decay_solid": gamma_s, + }, +} + + +def temp_z_decay(t, rho_w, rho_s, c_w, c_s, gamma_w, gamma_s, n, T0): + t = np.atleast_1d(t) + coeff = (-gamma_w * n - gamma_s * (1 - n) * rho_s) / ( + rho_w * c_w * n + rho_s * c_s * (1 - n) + ) + return coeff * t + T0 + + +def build_models(idx, test): + # Base MF6 GWF model type + ws = test.workspace + name = cases[idx] + gwename = "gwe-" + name + + decay_kwargs = parameters[name] + + print(f"Building MF6 model...{name}") + + sim = flopy.mf6.MFSimulation( + sim_name="heat", + sim_ws=ws, + exe_name="mf6", + version="mf6", + ) + tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(perlen, nstp, 1.0)]) + ims = flopy.mf6.ModflowIms( + sim, complexity="SIMPLE", inner_dvclose=0.001 + ) # T can not become negative in this model + + gwe = flopy.mf6.ModflowGwe( + sim, + modelname=gwename, + save_flows=True, + model_nam_file=f"{gwename}.nam", + ) + dis = flopy.mf6.ModflowGwedis( + gwe, nrow=nrow, ncol=ncol, nlay=nlay, delr=delr, delc=delc, top=top, botm=botm + ) + ic = flopy.mf6.ModflowGweic(gwe, strt=T0) + est = flopy.mf6.ModflowGweest( + gwe, + zero_order_decay_water=decay_kwargs["zero_order_decay_water"], + zero_order_decay_solid=decay_kwargs["zero_order_decay_solid"], + density_water=rho_w, + density_solid=rho_s, + heat_capacity_water=c_w, + heat_capacity_solid=c_s, + porosity=n, + decay_water=decay_kwargs["decay_water"], + decay_solid=decay_kwargs["decay_solid"], + ) + + oc = flopy.mf6.ModflowGweoc( + gwe, + budget_filerecord=f"{gwe.name}.bud", + temperature_filerecord=f"{gwe.name}.ucn", + printrecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")], + saverecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")], + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + + # read transport results from GWE model + name = cases[idx] + gwename = "gwe-" + name + + # Get the MF6 temperature output + sim = test.sims[0] + gwe = sim.get_model(gwename) + temp_ts = gwe.output.temperature().get_ts((0, 0, 0)) + t = temp_ts[:, 0] + + temp_analy_w = temp_z_decay( + t, rho_w, rho_s, c_w, c_s, gamma_w, 0, n, T0 + ) # aqueous decay only + temp_analy_s = temp_z_decay( + t, rho_w, rho_s, c_w, c_s, 0, gamma_s, n, T0 + ) # aqueous decay only + temp_analy_ws = temp_z_decay( + t, rho_w, rho_s, c_w, c_s, gamma_w, gamma_s, n, T0 + ) # aqueous + solid decay + + print("temperature evaluation: " + str(temp_analy_w)) + + if "aqe" in name: + assert temp_ts[-1, 1] == temp_analy_w[-1] + + if "sld" in name: + assert temp_ts[-1, 1] == temp_analy_s[-1] + + if "both" in name: + assert temp_ts[-1, 1] == temp_analy_ws[-1] + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 973dae2d8a9..4557d392d11 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -25,6 +25,7 @@ \item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed. \item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed. \item A profiling module is added for more enhanced performance diagnostics of the program. It can be activated through the PROFILE\_OPTION in the simulation name file. + \item Energy decay in the solid phase was added to the EST Package. This capability is described in the Supplemental Technical Information document, but no actual support was provided in the source code. Users can now specify distinct zeroth-order decay rates in either the aqueous or solid phases. \end{itemize} \underline{INTERNAL FLOW PACKAGES} diff --git a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn index 283900bd6b2..07a2ab1b100 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn @@ -9,12 +9,20 @@ longname save calculated flows to budget file description REPLACE save_flows {'{#1}': 'EST'} block options -name zero_order_decay +name zero_order_decay_water type keyword reader urword optional true -longname activate zero-order decay -description is a text keyword to indicate that zero-order decay will occur. Use of this keyword requires that DECAY is specified in the GRIDDATA block. +longname activate zero-order decay in aqueous phase +description is a text keyword to indicate that zero-order decay will occur in the aqueous phase. Use of this keyword requires that DECAY_WATER is specified in the GRIDDATA block. + +block options +name zero_order_decay_solid +type keyword +reader urword +optional true +longname activate zero-order decay in solid phase +description is a text keyword to indicate that zero-order decay will occur in the solid phase. Use of this keyword requires that DECAY_SOLID is specified in the GRIDDATA block. block options name density_water @@ -58,14 +66,24 @@ longname porosity description is the mobile domain porosity, defined as the mobile domain pore volume per mobile domain volume. The GWE model does not support the concept of an immobile domain in the context of heat transport. block griddata -name decay +name decay_water type double precision shape (nodes) reader readarray layered true optional true longname aqueous phase decay rate coefficient -description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay is energy per length cubed per time. Zero-order decay will have no effect on simulation results unless zero-order decay is specified in the options block. +description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO_ORDER_DECAY_WATER is specified in the options block. + +block griddata +name decay_solid +type double precision +shape (nodes) +reader readarray +layered true +optional true +longname solid phase decay rate coefficient +description is the rate coefficient for zero-order decay for the solid phase. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the solid phase is energy per length cubed per time. Zero-order decay in the solid phase will have no effect on simulation results unless ZERO_ORDER_DECAY_SOLID is specified in the options block. block griddata name heat_capacity_solid diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 4f6d3d27db2..cfe5784380b 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -13,7 +13,7 @@ module GweEstModule use KindModule, only: DP, I4B - use ConstantsModule, only: DONE, DZERO, DTWO, DHALF, LENBUDTXT, DEP3 + use ConstantsModule, only: DONE, IZERO, DZERO, DTWO, DHALF, LENBUDTXT, DEP3 use SimVariablesModule, only: errmsg, warnmsg use SimModule, only: store_error, count_errors, & store_warning @@ -27,9 +27,9 @@ module GweEstModule public :: GweEstType public :: est_cr ! - integer(I4B), parameter :: NBDITEMS = 2 + integer(I4B), parameter :: NBDITEMS = 3 character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt - data budtxt/' STORAGE-CELLBLK', ' DECAY-AQUEOUS'/ + data budtxt/' STORAGE-CELLBLK', ' DECAY-AQUEOUS', ' DECAY-SOLID'/ !> @ brief Energy storage and transfer !! @@ -47,10 +47,14 @@ module GweEstModule real(DP), dimension(:), pointer, contiguous :: ratesto => null() !< rate of energy storage ! ! -- decay - integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero) - real(DP), dimension(:), pointer, contiguous :: decay => null() !< first or zero order decay rate (aqueous) - real(DP), dimension(:), pointer, contiguous :: ratedcy => null() !< rate of decay - real(DP), dimension(:), pointer, contiguous :: decaylast => null() !< decay rate used for last iteration (needed for zero order decay) + integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero (aqueous and/or solid)) + integer(I4B), dimension(:), pointer :: idcytype => null() !< decay source (or sink) (0: not specified + real(DP), dimension(:), pointer, contiguous :: decay_water => null() !< first or zero order decay rate (aqueous) + real(DP), dimension(:), pointer, contiguous :: decay_solid => null() !< first or zero order decay rate (solid) + real(DP), dimension(:), pointer, contiguous :: ratedcyw => null() !< rate of decay in aqueous phase + real(DP), dimension(:), pointer, contiguous :: ratedcys => null() !< rate of decay in solid phase + real(DP), dimension(:), pointer, contiguous :: decaylastw => null() !< aqueous phase decay rate used for last iteration (needed for zero order decay) + real(DP), dimension(:), pointer, contiguous :: decaylasts => null() !< solid phase decay rate used for last iteration (needed for zero order decay) ! ! -- misc integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to model ibound @@ -64,9 +68,11 @@ module GweEstModule procedure :: est_fc procedure :: est_fc_sto procedure :: est_fc_dcy + procedure :: est_fc_dcy_solid procedure :: est_cq procedure :: est_cq_sto procedure :: est_cq_dcy + procedure :: est_cq_dcy_solid procedure :: est_bd procedure :: est_ot_flow procedure :: est_da @@ -156,7 +162,6 @@ end subroutine est_ar !< subroutine est_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, & rhs, kiter) - ! -- modules ! -- dummy class(GweEstType) :: this !< GweEstType object integer, intent(in) :: nodes !< number of nodes @@ -167,7 +172,6 @@ subroutine est_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, & real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model real(DP), intent(in), dimension(nodes) :: cnew !< temperature at end of this time step integer(I4B), intent(in) :: kiter !< solution outer iteration number - ! -- local ! ! -- storage contribution call this%est_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs) @@ -176,6 +180,8 @@ subroutine est_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, & if (this%idcy /= 0) then call this%est_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, & rhs, kiter) + call this%est_fc_dcy_solid(nodes, cold, nja, matrix_sln, idxglo, rhs, & + cnew, kiter) end if end subroutine est_fc @@ -268,18 +274,19 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & ! ! -- first order decay rate is a function of temperature, so add ! note: May want to remove first-order decay for temperature and support only zero-order ! to left hand side - hhcof = -this%decay(n) * vcell * swtpdt * this%porosity(n) & + hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac call matrix_sln%add_value_pos(idxglo(idiag), hhcof) - elseif (this%idcy == 2) then + elseif (this%idcy == 2 .and. this%idcytype(1) > 0) then ! ! -- Call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative temperatures ! Important note: still need to think through negative temps - decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), & - kiter, cold(n), cnew(n), delt) + decay_rate = get_zero_order_decay(this%decay_water(n), & + this%decaylastw(n), kiter, cold(n), & + cnew(n), delt) ! -- This term does get divided by eqnsclfac for fc purposes because it ! should start out being a rate of energy - this%decaylast(n) = decay_rate + this%decaylastw(n) = decay_rate rrhs = decay_rate * vcell * swtpdt * this%porosity(n) rhs(n) = rhs(n) + rrhs end if @@ -287,6 +294,62 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & end do end subroutine est_fc_dcy + !> @ brief Fill solid decay coefficient method for package + !! + !! Method to calculate and fill energy decay coefficients for the solid phase. + !< + subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & + rhs, cnew, kiter) + ! -- modules + use TdisModule, only: delt + ! -- dummy + class(GweEstType) :: this !< GwtMstType object + integer, intent(in) :: nodes !< number of nodes + real(DP), intent(in), dimension(nodes) :: cold !< temperature at end of last time step + integer(I4B), intent(in) :: nja !< number of GWE connections + class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix + integer(I4B), intent(in), dimension(nja) :: idxglo !< mapping vector for model (local) to solution (global) + real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model + real(DP), intent(in), dimension(nodes) :: cnew !< temperature at end of this time step + integer(I4B), intent(in) :: kiter !< solution outer iteration number + ! -- local + integer(I4B) :: n, idiag + real(DP) :: hhcof, rrhs + real(DP) :: vcell + real(DP) :: decay_rate + ! + ! -- loop through and calculate sorption contribution to hcof and rhs + do n = 1, this%dis%nodes + ! + ! -- skip if transport inactive + if (this%ibound(n) <= 0) cycle + ! + ! -- set variables + hhcof = DZERO + rrhs = DZERO + vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) + ! + ! -- add decay rate terms for the solid phase to the accumulators + idiag = this%dis%con%ia(n) + ! + ! -- add sorbed mass decay rate terms to accumulators + if (this%idcy == 1) then + ! -- first order decay rate not supported in GWE + elseif (this%idcy == 2 .and. this%idcytype(2) > 0) then + ! + ! -- Call function to get zero-order decay rate for the solid phase, + ! which may be changed from the user-specified rate to prevent + ! negative concentrations + decay_rate = get_zero_order_decay(this%decay_solid(n), & + this%decaylasts(n), & + kiter, cold(n), cnew(n), delt) + this%decaylasts(n) = decay_rate + rrhs = decay_rate * vcell * (1 - this%porosity(n)) * this%rhos(n) + rhs(n) = rhs(n) + rrhs + end if + end do + end subroutine est_fc_dcy_solid + !> @ brief Calculate flows for package !! !! Method to calculate flows for the package. @@ -306,7 +369,12 @@ subroutine est_cq(this, nodes, cnew, cold, flowja) ! ! -- decay if (this%idcy /= 0) then - call this%est_cq_dcy(nodes, cnew, cold, flowja) + if (this%idcytype(1) > 0) then + call this%est_cq_dcy(nodes, cnew, cold, flowja) + end if + if (this%idcytype(2) > 0) then + call this%est_cq_dcy_solid(nodes, cnew, cold, flowja) + end if end if end subroutine est_cq @@ -362,9 +430,9 @@ subroutine est_cq_sto(this, nodes, cnew, cold, flowja) end do end subroutine est_cq_sto - !> @ brief Calculate decay terms for package + !> @ brief Calculate decay terms for aqueous phase !! - !! Method to calculate decay terms for the package. + !! Method to calculate decay terms for the aqueous phase. !< subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this handles only decay in water; need to add zero-order (but not first-order?) decay in solid ! -- modules @@ -390,7 +458,7 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha do n = 1, nodes ! ! -- skip if transport inactive - this%ratedcy(n) = DZERO + this%ratedcyw(n) = DZERO if (this%ibound(n) <= 0) cycle ! ! -- calculate new and old water volumes @@ -402,21 +470,76 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha hhcof = DZERO rrhs = DZERO if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? - hhcof = -this%decay(n) * vcell * swtpdt * this%porosity(n) & + hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac - elseif (this%idcy == 2) then - decay_rate = get_zero_order_decay(this%decay(n), this%decaylast(n), & - 0, cold(n), cnew(n), delt) + elseif (this%idcy == 2 .and. this%idcytype(1) > IZERO) then ! zero order decay aqueous phase + decay_rate = get_zero_order_decay(this%decay_water(n), & + this%decaylastw(n), 0, cold(n), & + cnew(n), delt) rrhs = decay_rate * vcell * swtpdt * this%porosity(n) ! Important note: this term does NOT get multiplied by eqnsclfac for cq purposes because it should already be a rate of energy end if rate = hhcof * cnew(n) - rrhs - this%ratedcy(n) = rate + this%ratedcyw(n) = rate idiag = this%dis%con%ia(n) flowja(idiag) = flowja(idiag) + rate ! end do end subroutine est_cq_dcy + !> @ brief Calculate decay terms for aqueous phase + !! + !! Method to calculate decay terms for the aqueous phase. + !< + subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: this handles only decay in solid + ! -- modules + use TdisModule, only: delt + ! -- dummy + class(GweEstType) :: this !< GweEstType object + integer(I4B), intent(in) :: nodes !< number of nodes + real(DP), intent(in), dimension(nodes) :: cnew !< temperature at end of this time step + real(DP), intent(in), dimension(nodes) :: cold !< temperature at end of last time step + real(DP), dimension(:), contiguous, intent(inout) :: flowja !< flow between two connected control volumes + ! -- local + integer(I4B) :: n + integer(I4B) :: idiag + real(DP) :: rate + real(DP) :: hhcof, rrhs + real(DP) :: vcell + real(DP) :: decay_rate + ! + ! -- initialize + ! + ! -- Calculate decay change + do n = 1, nodes + ! + ! -- skip if transport inactive + this%ratedcys(n) = DZERO + if (this%ibound(n) <= 0) cycle + ! + ! -- calculate new and old water volumes + vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) + ! + ! -- calculate decay gains and losses + rate = DZERO + hhcof = DZERO + rrhs = DZERO + if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? + hhcof = -this%decay_solid(n) * vcell * (1 - this%porosity(n)) & + * this%eqnsclfac + elseif (this%idcy == 2 .and. this%idcytype(2) > IZERO) then ! zero order decay in the solid phase + decay_rate = get_zero_order_decay(this%decay_solid(n), & + this%decaylasts(n), 0, cold(n), & + cnew(n), delt) + rrhs = decay_rate * vcell * (1 - this%porosity(n)) * this%rhos(n) ! Important note: this term does NOT get multiplied by eqnsclfac for cq purposes because it should already be a rate of energy + end if + rate = hhcof * cnew(n) - rrhs + this%ratedcys(n) = rate + idiag = this%dis%con%ia(n) + flowja(idiag) = flowja(idiag) + rate + ! + end do + end subroutine est_cq_dcy_solid + !> @ brief Calculate budget terms for package !! !! Method to calculate budget terms for the package. @@ -440,9 +563,18 @@ subroutine est_bd(this, isuppress_output, model_budget) ! ! -- dcy if (this%idcy /= 0) then - call rate_accumulator(this%ratedcy, rin, rout) - call model_budget%addentry(rin, rout, delt, budtxt(2), & - isuppress_output, rowlabel=this%packName) + if (this%idcytype(1) > IZERO) then + ! -- aqueous phase + call rate_accumulator(this%ratedcyw, rin, rout) + call model_budget%addentry(rin, rout, delt, budtxt(2), & + isuppress_output, rowlabel=this%packName) + end if + if (this%idcytype(2) > IZERO) then + ! -- solid phase + call rate_accumulator(this%ratedcys, rin, rout) + call model_budget%addentry(rin, rout, delt, budtxt(3), & + isuppress_output, rowlabel=this%packName) + end if end if end subroutine est_bd @@ -457,7 +589,6 @@ subroutine est_ot_flow(this, icbcfl, icbcun) integer(I4B), intent(in) :: icbcun !< flag indication if cell-by-cell data should be saved ! -- local integer(I4B) :: ibinun - !character(len=16), dimension(2) :: aname integer(I4B) :: iprint, nvaluesp, nwidthp character(len=1) :: cdatafmp = ' ', editdesc = ' ' real(DP) :: dinact @@ -483,10 +614,20 @@ subroutine est_ot_flow(this, icbcfl, icbcun) nwidthp, editdesc, dinact) ! ! -- dcy - if (this%idcy /= 0) & - call this%dis%record_array(this%ratedcy, this%iout, iprint, -ibinun, & - budtxt(2), cdatafmp, nvaluesp, & - nwidthp, editdesc, dinact) + if (this%idcy /= 0) then + if (this%idcytype(1) > IZERO) then + ! -- aqueous phase + call this%dis%record_array(this%ratedcyw, this%iout, iprint, & + -ibinun, budtxt(2), cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + end if + if (this%idcytype(2) > IZERO) then + ! -- solid phase + call this%dis%record_array(this%ratedcys, this%iout, iprint, & + -ibinun, budtxt(3), cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + end if + end if end if end subroutine est_ot_flow @@ -505,9 +646,13 @@ subroutine est_da(this) call mem_deallocate(this%porosity) call mem_deallocate(this%ratesto) call mem_deallocate(this%idcy) - call mem_deallocate(this%decay) - call mem_deallocate(this%ratedcy) - call mem_deallocate(this%decaylast) + call mem_deallocate(this%idcytype) + call mem_deallocate(this%decay_water) + call mem_deallocate(this%decay_solid) + call mem_deallocate(this%ratedcyw) + call mem_deallocate(this%ratedcys) + call mem_deallocate(this%decaylastw) + call mem_deallocate(this%decaylasts) call mem_deallocate(this%cpw) call mem_deallocate(this%cps) call mem_deallocate(this%rhow) @@ -533,6 +678,7 @@ subroutine allocate_scalars(this) ! -- dummy class(GweEstType) :: this !< GweEstType object ! -- local + integer(I4B) :: n ! ! -- Allocate scalars in NumericalPackageType call this%NumericalPackageType%allocate_scalars() @@ -542,12 +688,16 @@ subroutine allocate_scalars(this) call mem_allocate(this%rhow, 'RHOW', this%memoryPath) call mem_allocate(this%latheatvap, 'LATHEATVAP', this%memoryPath) call mem_allocate(this%idcy, 'IDCY', this%memoryPath) + call mem_allocate(this%idcytype, 2, 'IDCYTYPE', this%memoryPath) ! ! -- Initialize this%cpw = DZERO this%rhow = DZERO this%latheatvap = DZERO this%idcy = 0 + do n = 1, 2 + this%idcytype(n) = IZERO + end do end subroutine allocate_scalars !> @ brief Allocate arrays for package @@ -573,13 +723,21 @@ subroutine allocate_arrays(this, nodes) ! ! -- dcy if (this%idcy == 0) then - call mem_allocate(this%ratedcy, 1, 'RATEDCY', this%memoryPath) - call mem_allocate(this%decay, 1, 'DECAY', this%memoryPath) - call mem_allocate(this%decaylast, 1, 'DECAYLAST', this%memoryPath) + call mem_allocate(this%ratedcyw, 1, 'RATEDCYW', this%memoryPath) + call mem_allocate(this%ratedcys, 1, 'RATEDCYS', this%memoryPath) + call mem_allocate(this%decay_water, 1, 'DECAY_WATER', this%memoryPath) + call mem_allocate(this%decay_solid, 1, 'DECAY_SOLID', this%memoryPath) + call mem_allocate(this%decaylastw, 1, 'DECAYLASTW', this%memoryPath) + call mem_allocate(this%decaylasts, 1, 'DECAYLAST', this%memoryPath) else - call mem_allocate(this%ratedcy, this%dis%nodes, 'RATEDCY', this%memoryPath) - call mem_allocate(this%decay, nodes, 'DECAY', this%memoryPath) - call mem_allocate(this%decaylast, nodes, 'DECAYLAST', this%memoryPath) + call mem_allocate(this%ratedcyw, this%dis%nodes, 'RATEDCYW', & + this%memoryPath) + call mem_allocate(this%ratedcys, this%dis%nodes, 'RATEDCYS', & + this%memoryPath) + call mem_allocate(this%decay_water, nodes, 'DECAY_WATER', this%memoryPath) + call mem_allocate(this%decay_solid, nodes, 'DECAY_SOLID', this%memoryPath) + call mem_allocate(this%decaylastw, nodes, 'DECAYLASTW', this%memoryPath) + call mem_allocate(this%decaylasts, nodes, 'DECAYLASTS', this%memoryPath) end if ! ! -- Initialize @@ -589,10 +747,13 @@ subroutine allocate_arrays(this, nodes) this%cps(n) = DZERO this%rhos(n) = DZERO end do - do n = 1, size(this%decay) - this%decay(n) = DZERO - this%ratedcy(n) = DZERO - this%decaylast(n) = DZERO + do n = 1, size(this%decay_water) + this%decay_water(n) = DZERO + this%decay_solid(n) = DZERO + this%ratedcyw(n) = DZERO + this%ratedcys(n) = DZERO + this%decaylastw(n) = DZERO + this%decaylasts(n) = DZERO end do end subroutine allocate_arrays @@ -616,7 +777,11 @@ subroutine read_options(this) character(len=*), parameter :: fmtidcy1 = & "(4x,'FIRST-ORDER DECAY IS ACTIVE. ')" character(len=*), parameter :: fmtidcy2 = & - "(4x,'ZERO-ORDER DECAY IS ACTIVE. ')" + "(4x,'ZERO-ORDER DECAY IN THE AQUEOUS "// & + &"PHASE IS ACTIVE. ')" + character(len=*), parameter :: fmtidcy3 = & + "(4x,'ZERO-ORDER DECAY IN THE SOLID "// & + &"PHASE IS ACTIVE. ')" ! ! -- get options block call this%parser%GetBlock('OPTIONS', isfound, ierr, blockRequired=.false., & @@ -636,9 +801,14 @@ subroutine read_options(this) case ('FIRST_ORDER_DECAY') this%idcy = 1 write (this%iout, fmtidcy1) - case ('ZERO_ORDER_DECAY') + case ('ZERO_ORDER_DECAY_WATER') this%idcy = 2 + this%idcytype(1) = 1 write (this%iout, fmtidcy2) + case ('ZERO_ORDER_DECAY_SOLID') + this%idcy = 2 + this%idcytype(2) = 2 + write (this%iout, fmtidcy3) case ('HEAT_CAPACITY_WATER') this%cpw = this%parser%GetDouble() if (this%cpw <= 0.0) then @@ -693,14 +863,15 @@ subroutine read_data(this) character(len=:), allocatable :: line integer(I4B) :: istart, istop, lloc, ierr logical :: isfound, endOfBlock - logical, dimension(4) :: lname - character(len=24), dimension(4) :: aname + logical, dimension(5) :: lname + character(len=24), dimension(5) :: aname ! -- formats ! -- data data aname(1)/' MOBILE DOMAIN POROSITY'/ - data aname(2)/' DECAY RATE'/ - data aname(3)/' HEAT CAPACITY OF SOLIDS'/ - data aname(4)/' DENSITY OF SOLIDS'/ + data aname(2)/'DECAY RATE AQUEOUS PHASE'/ + data aname(3)/' DECAY RATE SOLID PHASE'/ + data aname(4)/' HEAT CAPACITY OF SOLIDS'/ + data aname(5)/' DENSITY OF SOLIDS'/ ! ! -- initialize isfound = .false. @@ -722,24 +893,32 @@ subroutine read_data(this) this%parser%iuactive, this%porosity, & aname(1)) lname(1) = .true. - case ('DECAY') + case ('DECAY_WATER') if (this%idcy == 0) & - call mem_reallocate(this%decay, this%dis%nodes, 'DECAY', & + call mem_reallocate(this%decay_water, this%dis%nodes, 'DECAY_WATER', & trim(this%memoryPath)) call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & - this%parser%iuactive, this%decay, & + this%parser%iuactive, this%decay_water, & aname(2)) lname(2) = .true. - case ('HEAT_CAPACITY_SOLID') + case ('DECAY_SOLID') + if (this%idcy == 0) & + call mem_reallocate(this%decay_solid, this%dis%nodes, 'DECAY_SOLID', & + trim(this%memoryPath)) call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & - this%parser%iuactive, this%cps, & + this%parser%iuactive, this%decay_solid, & aname(3)) lname(3) = .true. - case ('DENSITY_SOLID') + case ('HEAT_CAPACITY_SOLID') call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & - this%parser%iuactive, this%rhos, & + this%parser%iuactive, this%cps, & aname(4)) lname(4) = .true. + case ('DENSITY_SOLID') + call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & + this%parser%iuactive, this%rhos, & + aname(5)) + lname(5) = .true. case default write (errmsg, '(a,a)') 'Unknown griddata tag: ', trim(keyword) call store_error(errmsg) @@ -758,28 +937,38 @@ subroutine read_data(this) write (errmsg, '(a)') 'Porosity not specified in griddata block.' call store_error(errmsg) end if - if (.not. lname(3)) then + if (.not. lname(4)) then write (errmsg, '(a)') 'HEAT_CAPACITY_SOLID not specified in griddata block.' call store_error(errmsg) end if - if (.not. lname(4)) then + if (.not. lname(5)) then write (errmsg, '(a)') 'DENSITY_SOLID not specified in griddata block.' call store_error(errmsg) end if ! ! -- Check for required decay/production rate coefficients if (this%idcy > 0) then - if (.not. lname(2)) then - write (errmsg, '(a)') 'First or zero order decay is & - &active but the first rate coefficient is not specified. Decay & - &must be specified in griddata block.' + if (.not. lname(2) .and. .not. lname(3)) then + write (errmsg, '(a)') 'Zero order decay in either the aqueous & + &or solid phase is active but zero-order rate coefficients, & + &either in the aqueous or solid phase, is not specified. & + &Either DECAY_WATER or DECAY_SOLID must be specified in the & + &griddata block.' call store_error(errmsg) end if else if (lname(2)) then - write (warnmsg, '(a)') 'First or zero orer decay & - &is not active but decay was specified. Decay will & - &have no affect on simulation results.' + write (warnmsg, '(a)') 'Zero order decay in the aqueous phase has & + ¬ been activated but DECAY_WATER has been specified. Zero & + &order decay in the aqueous phase will have no affect on & + &simulation results.' + call store_warning(warnmsg) + write (this%iout, '(1x,a)') 'WARNING. '//warnmsg + else if (lname(3)) then + write (warnmsg, '(a)') 'Zero order decay in the solid phase has not & + &been activated but DECAY_SOLID has been specified. Zero order & + &decay in the solid phase will have no affect on simulation & + &results.' call store_warning(warnmsg) write (this%iout, '(1x,a)') 'WARNING. '//warnmsg end if diff --git a/src/Model/GroundWaterTransport/gwt-mst.f90 b/src/Model/GroundWaterTransport/gwt-mst.f90 index 91c0212ee9a..2f95d72d84c 100644 --- a/src/Model/GroundWaterTransport/gwt-mst.f90 +++ b/src/Model/GroundWaterTransport/gwt-mst.f90 @@ -441,7 +441,6 @@ end subroutine mst_srb_term !> @ brief Fill sorption-decay coefficient method for package !! !! Method to calculate and fill sorption-decay coefficients for the package. - !! !< subroutine mst_fc_dcy_srb(this, nodes, cold, nja, matrix_sln, idxglo, & rhs, cnew, kiter) From 8dfc1c34e07776c9c8b4312acaaeca2648e7f0ca Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 20 Jan 2025 19:33:50 -0800 Subject: [PATCH 02/18] different implementation of new variable idcytype --- src/Model/GroundWaterEnergy/gwe-est.f90 | 50 ++++++++++++++----------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index cfe5784380b..6582600a98f 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -48,7 +48,7 @@ module GweEstModule ! ! -- decay integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero (aqueous and/or solid)) - integer(I4B), dimension(:), pointer :: idcytype => null() !< decay source (or sink) (0: not specified + integer(I4B), pointer :: idcytype => null() !< decay source (or sink) (1: aqueous only, 2: solid only, 3: both phases real(DP), dimension(:), pointer, contiguous :: decay_water => null() !< first or zero order decay rate (aqueous) real(DP), dimension(:), pointer, contiguous :: decay_solid => null() !< first or zero order decay rate (solid) real(DP), dimension(:), pointer, contiguous :: ratedcyw => null() !< rate of decay in aqueous phase @@ -277,7 +277,8 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac call matrix_sln%add_value_pos(idxglo(idiag), hhcof) - elseif (this%idcy == 2 .and. this%idcytype(1) > 0) then + elseif (this%idcy == 2 .and. (this%idcytype == 1 .or. & + this%idcytype == 3)) then ! ! -- Call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative temperatures ! Important note: still need to think through negative temps @@ -335,7 +336,8 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & ! -- add sorbed mass decay rate terms to accumulators if (this%idcy == 1) then ! -- first order decay rate not supported in GWE - elseif (this%idcy == 2 .and. this%idcytype(2) > 0) then + elseif (this%idcy == 2 .and. (this%idcytype == 2 .or. & + this%idcytype == 3)) then ! ! -- Call function to get zero-order decay rate for the solid phase, ! which may be changed from the user-specified rate to prevent @@ -369,10 +371,10 @@ subroutine est_cq(this, nodes, cnew, cold, flowja) ! ! -- decay if (this%idcy /= 0) then - if (this%idcytype(1) > 0) then + if (this%idcytype == 1 .or. this%idcytype == 3) then call this%est_cq_dcy(nodes, cnew, cold, flowja) end if - if (this%idcytype(2) > 0) then + if (this%idcytype == 2 .or. this%idcytype == 3) then call this%est_cq_dcy_solid(nodes, cnew, cold, flowja) end if end if @@ -452,8 +454,6 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha real(DP) :: vcell real(DP) :: decay_rate ! - ! -- initialize - ! ! -- Calculate decay change do n = 1, nodes ! @@ -472,7 +472,8 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcytype(1) > IZERO) then ! zero order decay aqueous phase + elseif (this%idcy == 2 .and. this%idcytype == 1 .or. & + this%idcytype == 3) then ! zero order decay aqueous phase decay_rate = get_zero_order_decay(this%decay_water(n), & this%decaylastw(n), 0, cold(n), & cnew(n), delt) @@ -526,7 +527,8 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: t if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_solid(n) * vcell * (1 - this%porosity(n)) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcytype(2) > IZERO) then ! zero order decay in the solid phase + elseif (this%idcy == 2 .and. this%idcytype == 2 .or. & + this%idcytype == 3) then ! zero order decay in the solid phase decay_rate = get_zero_order_decay(this%decay_solid(n), & this%decaylasts(n), 0, cold(n), & cnew(n), delt) @@ -563,13 +565,13 @@ subroutine est_bd(this, isuppress_output, model_budget) ! ! -- dcy if (this%idcy /= 0) then - if (this%idcytype(1) > IZERO) then + if (this%idcytype == 1 .or. this%idcytype == 3) then ! -- aqueous phase call rate_accumulator(this%ratedcyw, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(2), & isuppress_output, rowlabel=this%packName) end if - if (this%idcytype(2) > IZERO) then + if (this%idcytype == 2 .or. this%idcytype == 3) then ! -- solid phase call rate_accumulator(this%ratedcys, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(3), & @@ -615,13 +617,13 @@ subroutine est_ot_flow(this, icbcfl, icbcun) ! ! -- dcy if (this%idcy /= 0) then - if (this%idcytype(1) > IZERO) then + if (this%idcytype == 1 .or. this%idcytype == 3) then ! -- aqueous phase call this%dis%record_array(this%ratedcyw, this%iout, iprint, & -ibinun, budtxt(2), cdatafmp, nvaluesp, & nwidthp, editdesc, dinact) end if - if (this%idcytype(2) > IZERO) then + if (this%idcytype == 2 .or. this%idcytype == 3) then ! -- solid phase call this%dis%record_array(this%ratedcys, this%iout, iprint, & -ibinun, budtxt(3), cdatafmp, nvaluesp, & @@ -677,8 +679,6 @@ subroutine allocate_scalars(this) use MemoryManagerModule, only: mem_allocate, mem_setptr ! -- dummy class(GweEstType) :: this !< GweEstType object - ! -- local - integer(I4B) :: n ! ! -- Allocate scalars in NumericalPackageType call this%NumericalPackageType%allocate_scalars() @@ -688,16 +688,14 @@ subroutine allocate_scalars(this) call mem_allocate(this%rhow, 'RHOW', this%memoryPath) call mem_allocate(this%latheatvap, 'LATHEATVAP', this%memoryPath) call mem_allocate(this%idcy, 'IDCY', this%memoryPath) - call mem_allocate(this%idcytype, 2, 'IDCYTYPE', this%memoryPath) + call mem_allocate(this%idcytype, 'IDCYTYPE', this%memoryPath) ! ! -- Initialize this%cpw = DZERO this%rhow = DZERO this%latheatvap = DZERO - this%idcy = 0 - do n = 1, 2 - this%idcytype(n) = IZERO - end do + this%idcy = IZERO + this%idcytype = IZERO end subroutine allocate_scalars !> @ brief Allocate arrays for package @@ -803,11 +801,19 @@ subroutine read_options(this) write (this%iout, fmtidcy1) case ('ZERO_ORDER_DECAY_WATER') this%idcy = 2 - this%idcytype(1) = 1 + if (this%idcytype > IZERO) then + this%idcytype = 3 + else + this%idcytype = 1 + end if write (this%iout, fmtidcy2) case ('ZERO_ORDER_DECAY_SOLID') this%idcy = 2 - this%idcytype(2) = 2 + if (this%idcytype > IZERO) then + this%idcytype = 3 + else + this%idcytype = 2 + end if write (this%iout, fmtidcy3) case ('HEAT_CAPACITY_WATER') this%cpw = this%parser%GetDouble() From dee2812fe801851d56687982831c0bf75a08b23d Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 20 Jan 2025 19:36:51 -0800 Subject: [PATCH 03/18] rename variable to something more indicative of what it is doing --- src/Model/GroundWaterEnergy/gwe-est.f90 | 48 ++++++++++++------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 6582600a98f..c18421e8bf0 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -48,7 +48,7 @@ module GweEstModule ! ! -- decay integer(I4B), pointer :: idcy => null() !< order of decay rate (0:none, 1:first, 2:zero (aqueous and/or solid)) - integer(I4B), pointer :: idcytype => null() !< decay source (or sink) (1: aqueous only, 2: solid only, 3: both phases + integer(I4B), pointer :: idcysrc => null() !< decay source (or sink) (1: aqueous only, 2: solid only, 3: both phases real(DP), dimension(:), pointer, contiguous :: decay_water => null() !< first or zero order decay rate (aqueous) real(DP), dimension(:), pointer, contiguous :: decay_solid => null() !< first or zero order decay rate (solid) real(DP), dimension(:), pointer, contiguous :: ratedcyw => null() !< rate of decay in aqueous phase @@ -277,8 +277,8 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac call matrix_sln%add_value_pos(idxglo(idiag), hhcof) - elseif (this%idcy == 2 .and. (this%idcytype == 1 .or. & - this%idcytype == 3)) then + elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & + this%idcysrc == 3)) then ! ! -- Call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative temperatures ! Important note: still need to think through negative temps @@ -336,8 +336,8 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & ! -- add sorbed mass decay rate terms to accumulators if (this%idcy == 1) then ! -- first order decay rate not supported in GWE - elseif (this%idcy == 2 .and. (this%idcytype == 2 .or. & - this%idcytype == 3)) then + elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. & + this%idcysrc == 3)) then ! ! -- Call function to get zero-order decay rate for the solid phase, ! which may be changed from the user-specified rate to prevent @@ -371,10 +371,10 @@ subroutine est_cq(this, nodes, cnew, cold, flowja) ! ! -- decay if (this%idcy /= 0) then - if (this%idcytype == 1 .or. this%idcytype == 3) then + if (this%idcysrc == 1 .or. this%idcysrc == 3) then call this%est_cq_dcy(nodes, cnew, cold, flowja) end if - if (this%idcytype == 2 .or. this%idcytype == 3) then + if (this%idcysrc == 2 .or. this%idcysrc == 3) then call this%est_cq_dcy_solid(nodes, cnew, cold, flowja) end if end if @@ -472,8 +472,8 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcytype == 1 .or. & - this%idcytype == 3) then ! zero order decay aqueous phase + elseif (this%idcy == 2 .and. this%idcysrc == 1 .or. & + this%idcysrc == 3) then ! zero order decay aqueous phase decay_rate = get_zero_order_decay(this%decay_water(n), & this%decaylastw(n), 0, cold(n), & cnew(n), delt) @@ -527,8 +527,8 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: t if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_solid(n) * vcell * (1 - this%porosity(n)) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcytype == 2 .or. & - this%idcytype == 3) then ! zero order decay in the solid phase + elseif (this%idcy == 2 .and. this%idcysrc == 2 .or. & + this%idcysrc == 3) then ! zero order decay in the solid phase decay_rate = get_zero_order_decay(this%decay_solid(n), & this%decaylasts(n), 0, cold(n), & cnew(n), delt) @@ -565,13 +565,13 @@ subroutine est_bd(this, isuppress_output, model_budget) ! ! -- dcy if (this%idcy /= 0) then - if (this%idcytype == 1 .or. this%idcytype == 3) then + if (this%idcysrc == 1 .or. this%idcysrc == 3) then ! -- aqueous phase call rate_accumulator(this%ratedcyw, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(2), & isuppress_output, rowlabel=this%packName) end if - if (this%idcytype == 2 .or. this%idcytype == 3) then + if (this%idcysrc == 2 .or. this%idcysrc == 3) then ! -- solid phase call rate_accumulator(this%ratedcys, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(3), & @@ -617,13 +617,13 @@ subroutine est_ot_flow(this, icbcfl, icbcun) ! ! -- dcy if (this%idcy /= 0) then - if (this%idcytype == 1 .or. this%idcytype == 3) then + if (this%idcysrc == 1 .or. this%idcysrc == 3) then ! -- aqueous phase call this%dis%record_array(this%ratedcyw, this%iout, iprint, & -ibinun, budtxt(2), cdatafmp, nvaluesp, & nwidthp, editdesc, dinact) end if - if (this%idcytype == 2 .or. this%idcytype == 3) then + if (this%idcysrc == 2 .or. this%idcysrc == 3) then ! -- solid phase call this%dis%record_array(this%ratedcys, this%iout, iprint, & -ibinun, budtxt(3), cdatafmp, nvaluesp, & @@ -648,7 +648,7 @@ subroutine est_da(this) call mem_deallocate(this%porosity) call mem_deallocate(this%ratesto) call mem_deallocate(this%idcy) - call mem_deallocate(this%idcytype) + call mem_deallocate(this%idcysrc) call mem_deallocate(this%decay_water) call mem_deallocate(this%decay_solid) call mem_deallocate(this%ratedcyw) @@ -688,14 +688,14 @@ subroutine allocate_scalars(this) call mem_allocate(this%rhow, 'RHOW', this%memoryPath) call mem_allocate(this%latheatvap, 'LATHEATVAP', this%memoryPath) call mem_allocate(this%idcy, 'IDCY', this%memoryPath) - call mem_allocate(this%idcytype, 'IDCYTYPE', this%memoryPath) + call mem_allocate(this%idcysrc, 'IDCYSRC', this%memoryPath) ! ! -- Initialize this%cpw = DZERO this%rhow = DZERO this%latheatvap = DZERO this%idcy = IZERO - this%idcytype = IZERO + this%idcysrc = IZERO end subroutine allocate_scalars !> @ brief Allocate arrays for package @@ -801,18 +801,18 @@ subroutine read_options(this) write (this%iout, fmtidcy1) case ('ZERO_ORDER_DECAY_WATER') this%idcy = 2 - if (this%idcytype > IZERO) then - this%idcytype = 3 + if (this%idcysrc > IZERO) then + this%idcysrc = 3 else - this%idcytype = 1 + this%idcysrc = 1 end if write (this%iout, fmtidcy2) case ('ZERO_ORDER_DECAY_SOLID') this%idcy = 2 - if (this%idcytype > IZERO) then - this%idcytype = 3 + if (this%idcysrc > IZERO) then + this%idcysrc = 3 else - this%idcytype = 2 + this%idcysrc = 2 end if write (this%iout, fmtidcy3) case ('HEAT_CAPACITY_WATER') From 378ed1575e6cad5de681fb469b46b45dea85902b Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 20 Jan 2025 19:57:56 -0800 Subject: [PATCH 04/18] Set a tolerance for value comparisons to 1e-10 --- autotest/test_gwe_decay01.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/autotest/test_gwe_decay01.py b/autotest/test_gwe_decay01.py index 7999b6b1e02..c2e324d1e5e 100644 --- a/autotest/test_gwe_decay01.py +++ b/autotest/test_gwe_decay01.py @@ -128,6 +128,13 @@ def build_models(idx, test): def check_output(idx, test): print("evaluating results...") + msg = ( + "Differences detected between the simulated results for zeroth-order " + "energy decay and the expected solution for decay specified in " + ) + msg0 = msg + "the aqueous phase." + msg1 = msg + "the solid phase." + msg2 = msg + "both the aqueous and solid phases." # read transport results from GWE model name = cases[idx] @@ -152,13 +159,13 @@ def check_output(idx, test): print("temperature evaluation: " + str(temp_analy_w)) if "aqe" in name: - assert temp_ts[-1, 1] == temp_analy_w[-1] + assert np.isclose(temp_ts[-1, 1], temp_analy_w[-1], atol=1e-10), msg0 if "sld" in name: - assert temp_ts[-1, 1] == temp_analy_s[-1] + assert np.isclose(temp_ts[-1, 1], temp_analy_s[-1], atol=1e-10), msg1 if "both" in name: - assert temp_ts[-1, 1] == temp_analy_ws[-1] + assert np.isclose(temp_ts[-1, 1], temp_analy_ws[-1], atol=1e-10), msg2 # - No need to change any code below From e9c1cbdc0f690bc3e195d8cd9fb7ee67c666fd33 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Mon, 20 Jan 2025 21:22:39 -0800 Subject: [PATCH 05/18] add latex escape character before underscores --- doc/mf6io/mf6ivar/dfn/gwe-est.dfn | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn index 07a2ab1b100..6172f361c66 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn @@ -14,7 +14,7 @@ type keyword reader urword optional true longname activate zero-order decay in aqueous phase -description is a text keyword to indicate that zero-order decay will occur in the aqueous phase. Use of this keyword requires that DECAY_WATER is specified in the GRIDDATA block. +description is a text keyword to indicate that zero-order decay will occur in the aqueous phase. Use of this keyword requires that DECAY\_WATER is specified in the GRIDDATA block. block options name zero_order_decay_solid @@ -22,7 +22,7 @@ type keyword reader urword optional true longname activate zero-order decay in solid phase -description is a text keyword to indicate that zero-order decay will occur in the solid phase. Use of this keyword requires that DECAY_SOLID is specified in the GRIDDATA block. +description is a text keyword to indicate that zero-order decay will occur in the solid phase. Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block. block options name density_water @@ -73,7 +73,7 @@ reader readarray layered true optional true longname aqueous phase decay rate coefficient -description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO_ORDER_DECAY_WATER is specified in the options block. +description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block. block griddata name decay_solid @@ -83,7 +83,7 @@ reader readarray layered true optional true longname solid phase decay rate coefficient -description is the rate coefficient for zero-order decay for the solid phase. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the solid phase is energy per length cubed per time. Zero-order decay in the solid phase will have no effect on simulation results unless ZERO_ORDER_DECAY_SOLID is specified in the options block. +description is the rate coefficient for zero-order decay for the solid phase. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the solid phase is energy per length cubed per time. Zero-order decay in the solid phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_SOLID is specified in the options block. block griddata name heat_capacity_solid From ae887bd6b01524d4a7f06233fa26b65ef9ff2910 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 14:58:06 -0800 Subject: [PATCH 06/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930711610 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index c18421e8bf0..2eccf0ad5f8 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -472,8 +472,8 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcysrc == 1 .or. & - this%idcysrc == 3) then ! zero order decay aqueous phase + elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & + this%idcysrc == 3)) then ! zero order decay aqueous phase decay_rate = get_zero_order_decay(this%decay_water(n), & this%decaylastw(n), 0, cold(n), & cnew(n), delt) @@ -527,8 +527,8 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: t if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? hhcof = -this%decay_solid(n) * vcell * (1 - this%porosity(n)) & * this%eqnsclfac - elseif (this%idcy == 2 .and. this%idcysrc == 2 .or. & - this%idcysrc == 3) then ! zero order decay in the solid phase + elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. & + this%idcysrc == 3)) then ! zero order decay in the solid phase decay_rate = get_zero_order_decay(this%decay_solid(n), & this%decaylasts(n), 0, cold(n), & cnew(n), delt) From a853d934d77dcdcdef650571dc2201ac884a9bda Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 15:19:38 -0800 Subject: [PATCH 07/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930697941 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 85 ++++++------------------- 1 file changed, 18 insertions(+), 67 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 2eccf0ad5f8..f1b25065419 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -282,9 +282,7 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & ! ! -- Call function to get zero-order decay rate, which may be changed ! from the user-specified rate to prevent negative temperatures ! Important note: still need to think through negative temps - decay_rate = get_zero_order_decay(this%decay_water(n), & - this%decaylastw(n), kiter, cold(n), & - cnew(n), delt) + decay_rate = this%decay_water(n) ! -- This term does get divided by eqnsclfac for fc purposes because it ! should start out being a rate of energy this%decaylastw(n) = decay_rate @@ -339,12 +337,12 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. & this%idcysrc == 3)) then ! - ! -- Call function to get zero-order decay rate for the solid phase, - ! which may be changed from the user-specified rate to prevent - ! negative concentrations - decay_rate = get_zero_order_decay(this%decay_solid(n), & - this%decaylasts(n), & - kiter, cold(n), cnew(n), delt) + ! -- negative temps are currently not checked for or prevented since a + ! user can define a temperature scale of their own choosing. if + ! negative temps result from the specified zero-order decay value, + ! it is up to the user to decide if the calculated temperatures are + ! acceptable + decay_rate = this%decay_solid(n) this%decaylasts(n) = decay_rate rrhs = decay_rate * vcell * (1 - this%porosity(n)) * this%rhos(n) rhs(n) = rhs(n) + rrhs @@ -474,10 +472,10 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha * this%eqnsclfac elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & this%idcysrc == 3)) then ! zero order decay aqueous phase - decay_rate = get_zero_order_decay(this%decay_water(n), & - this%decaylastw(n), 0, cold(n), & - cnew(n), delt) - rrhs = decay_rate * vcell * swtpdt * this%porosity(n) ! Important note: this term does NOT get multiplied by eqnsclfac for cq purposes because it should already be a rate of energy + decay_rate = this%decay_water(n) + ! -- this term does NOT get multiplied by eqnsclfac for cq purposes + ! because it should already be a rate of energy + rrhs = decay_rate * vcell * swtpdt * this%porosity(n) end if rate = hhcof * cnew(n) - rrhs this%ratedcyw(n) = rate @@ -524,15 +522,13 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: t rate = DZERO hhcof = DZERO rrhs = DZERO - if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? - hhcof = -this%decay_solid(n) * vcell * (1 - this%porosity(n)) & - * this%eqnsclfac - elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. & - this%idcysrc == 3)) then ! zero order decay in the solid phase - decay_rate = get_zero_order_decay(this%decay_solid(n), & - this%decaylasts(n), 0, cold(n), & - cnew(n), delt) - rrhs = decay_rate * vcell * (1 - this%porosity(n)) * this%rhos(n) ! Important note: this term does NOT get multiplied by eqnsclfac for cq purposes because it should already be a rate of energy + ! -- first-order decay (idcy=1) is not supported for temperature modeling + if (this%idcy == 2 .and. & + (this%idcysrc == 2 .or. this%idcysrc == 3)) then ! zero order decay in the solid phase + decay_rate = this%decay_solid(n) + ! -- this term does NOT get multiplied by eqnsclfac for cq purposes + ! because it should already be a rate of energy + rrhs = decay_rate * vcell * (1 - this%porosity(n)) * this%rhos(n) end if rate = hhcof * cnew(n) - rrhs this%ratedcys(n) = rate @@ -986,49 +982,4 @@ subroutine read_data(this) end if end subroutine read_data - !> @ brief Calculate zero-order decay rate and constrain if necessary - !! - !! Function to calculate the zero-order decay rate from the user specified - !! decay rate. If the decay rate is positive, then the decay rate must - !! be constrained so that more energy is not removed than is available. - !! Without this constraint, negative temperatures could result from - !! zero-order decay (no freezing). - !< - function get_zero_order_decay(decay_rate_usr, decay_rate_last, kiter, & - cold, cnew, delt) result(decay_rate) - ! -- dummy - real(DP), intent(in) :: decay_rate_usr !< user-entered decay rate - real(DP), intent(in) :: decay_rate_last !< decay rate used for last iteration - integer(I4B), intent(in) :: kiter !< Picard iteration counter - real(DP), intent(in) :: cold !< temperature at end of last time step - real(DP), intent(in) :: cnew !< temperature at end of this time step - real(DP), intent(in) :: delt !< length of time step - ! -- return - real(DP) :: decay_rate !< returned value for decay rate - ! - ! -- Return user rate if production, otherwise constrain, if necessary - if (decay_rate_usr < DZERO) then - ! - ! -- Production, no need to limit rate - decay_rate = decay_rate_usr - else - ! - ! -- Need to ensure decay does not result in negative - ! temperature, so reduce the rate if it would result in - ! removing more energy than is in the cell. ! kluge note: think through - if (kiter == 1) then - decay_rate = min(decay_rate_usr, cold / delt) ! kluge note: actually want to use rhow*cpw*cold and rhow*cpw*cnew for rates here and below - else - decay_rate = decay_rate_last - if (cnew < DZERO) then - decay_rate = decay_rate_last + cnew / delt - else if (cnew > cold) then - decay_rate = decay_rate_last + cnew / delt - end if - decay_rate = min(decay_rate_usr, decay_rate) - end if - decay_rate = max(decay_rate, DZERO) - end if - end function get_zero_order_decay - end module GweEstModule From 1f517175f417984187171e13af0c24ff64eab28d Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 15:34:37 -0800 Subject: [PATCH 08/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930666834 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index f1b25065419..16169af09a3 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -280,8 +280,6 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & this%idcysrc == 3)) then ! - ! -- Call function to get zero-order decay rate, which may be changed - ! from the user-specified rate to prevent negative temperatures ! Important note: still need to think through negative temps decay_rate = this%decay_water(n) ! -- This term does get divided by eqnsclfac for fc purposes because it ! should start out being a rate of energy @@ -952,16 +950,15 @@ subroutine read_data(this) if (this%idcy > 0) then if (.not. lname(2) .and. .not. lname(3)) then write (errmsg, '(a)') 'Zero order decay in either the aqueous & - &or solid phase is active but zero-order rate coefficients, & - &either in the aqueous or solid phase, is not specified. & - &Either DECAY_WATER or DECAY_SOLID must be specified in the & - &griddata block.' + &or solid phase is active but the corresponding zero-order & + &rate coefficient is not specified. Either DECAY_WATER or & + &DECAY_SOLID must be specified in the griddata block.' call store_error(errmsg) end if else if (lname(2)) then write (warnmsg, '(a)') 'Zero order decay in the aqueous phase has & - ¬ been activated but DECAY_WATER has been specified. Zero & + ¬ been activated but DECAY_WATER has been specified. Zero & &order decay in the aqueous phase will have no affect on & &simulation results.' call store_warning(warnmsg) From 8e3021e505a7c782326afa91e3fc0e50d28fbaed Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 15:40:09 -0800 Subject: [PATCH 09/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930616915 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 16169af09a3..e073de72e3c 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -239,8 +239,6 @@ end subroutine est_fc_sto !< subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & idxglo, rhs, kiter) - ! -- modules - use TdisModule, only: delt ! -- dummy class(GweEstType) :: this !< GweEstType object integer, intent(in) :: nodes !< number of nodes @@ -297,8 +295,6 @@ end subroutine est_fc_dcy !< subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & rhs, cnew, kiter) - ! -- modules - use TdisModule, only: delt ! -- dummy class(GweEstType) :: this !< GwtMstType object integer, intent(in) :: nodes !< number of nodes @@ -432,9 +428,7 @@ end subroutine est_cq_sto !! !! Method to calculate decay terms for the aqueous phase. !< - subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this handles only decay in water; need to add zero-order (but not first-order?) decay in solid - ! -- modules - use TdisModule, only: delt + subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! -- dummy class(GweEstType) :: this !< GweEstType object integer(I4B), intent(in) :: nodes !< number of nodes @@ -483,13 +477,11 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) ! Important note: this ha end do end subroutine est_cq_dcy - !> @ brief Calculate decay terms for aqueous phase + !> @ brief Calculate decay terms for solid phase !! - !! Method to calculate decay terms for the aqueous phase. + !! Method to calculate decay terms for the solid phase. !< - subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! Important note: this handles only decay in solid - ! -- modules - use TdisModule, only: delt + subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) ! -- dummy class(GweEstType) :: this !< GweEstType object integer(I4B), intent(in) :: nodes !< number of nodes From 6a93a3e7ad12932e4fd483f4249dca24b798dcae Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 15:54:10 -0800 Subject: [PATCH 10/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930623176 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 40 +++++++------------------ 1 file changed, 11 insertions(+), 29 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index e073de72e3c..279e0e61ac1 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -250,7 +250,7 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & real(DP), intent(inout), dimension(nodes) :: rhs !< right-hand side vector for model integer(I4B), intent(in) :: kiter !< solution outer iteration number ! -- local - integer(I4B) :: n, idiag + integer(I4B) :: n real(DP) :: hhcof, rrhs real(DP) :: swtpdt real(DP) :: vcell @@ -266,17 +266,9 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) swtpdt = this%fmi%gwfsat(n) ! - ! -- add decay rate terms to accumulators - idiag = this%dis%con%ia(n) - if (this%idcy == 1) then - ! - ! -- first order decay rate is a function of temperature, so add ! note: May want to remove first-order decay for temperature and support only zero-order - ! to left hand side - hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & - * this%eqnsclfac - call matrix_sln%add_value_pos(idxglo(idiag), hhcof) - elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & - this%idcysrc == 3)) then + ! -- add zero-order decay rate terms to accumulators + if (this%idcy == 2 .and. (this%idcysrc == 1 .or. & + this%idcysrc == 3)) then ! decay_rate = this%decay_water(n) ! -- This term does get divided by eqnsclfac for fc purposes because it @@ -306,7 +298,7 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & real(DP), intent(in), dimension(nodes) :: cnew !< temperature at end of this time step integer(I4B), intent(in) :: kiter !< solution outer iteration number ! -- local - integer(I4B) :: n, idiag + integer(I4B) :: n real(DP) :: hhcof, rrhs real(DP) :: vcell real(DP) :: decay_rate @@ -322,14 +314,9 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & rrhs = DZERO vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) ! - ! -- add decay rate terms for the solid phase to the accumulators - idiag = this%dis%con%ia(n) - ! - ! -- add sorbed mass decay rate terms to accumulators - if (this%idcy == 1) then - ! -- first order decay rate not supported in GWE - elseif (this%idcy == 2 .and. (this%idcysrc == 2 .or. & - this%idcysrc == 3)) then + ! -- account for zero-order decay rate terms in rhs + if (this%idcy == 2 .and. (this%idcysrc == 2 .or. & + this%idcysrc == 3)) then ! ! -- negative temps are currently not checked for or prevented since a ! user can define a temperature scale of their own choosing. if @@ -349,7 +336,6 @@ end subroutine est_fc_dcy_solid !! Method to calculate flows for the package. !< subroutine est_cq(this, nodes, cnew, cold, flowja) - ! -- modules ! -- dummy class(GweEstType) :: this !< GweEstType object integer(I4B), intent(in) :: nodes !< number of nodes @@ -459,11 +445,9 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) rate = DZERO hhcof = DZERO rrhs = DZERO - if (this%idcy == 1) then ! Important note: do we need/want first-order decay for temperature??? - hhcof = -this%decay_water(n) * vcell * swtpdt * this%porosity(n) & - * this%eqnsclfac - elseif (this%idcy == 2 .and. (this%idcysrc == 1 .or. & - this%idcysrc == 3)) then ! zero order decay aqueous phase + ! -- zero order decay aqueous phase + if (this%idcy == 2 .and. & + (this%idcysrc == 1 .or. this%idcysrc == 3)) then decay_rate = this%decay_water(n) ! -- this term does NOT get multiplied by eqnsclfac for cq purposes ! because it should already be a rate of energy @@ -496,8 +480,6 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) real(DP) :: vcell real(DP) :: decay_rate ! - ! -- initialize - ! ! -- Calculate decay change do n = 1, nodes ! From bf91985af76227e1bfecdcf91ce4bc464e05003c Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 16:04:21 -0800 Subject: [PATCH 11/18] Fix unused variable error --- src/Model/GroundWaterEnergy/gwe-est.f90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 279e0e61ac1..93b3b0a820f 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -251,7 +251,7 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & integer(I4B), intent(in) :: kiter !< solution outer iteration number ! -- local integer(I4B) :: n - real(DP) :: hhcof, rrhs + real(DP) :: rrhs real(DP) :: swtpdt real(DP) :: vcell real(DP) :: decay_rate @@ -299,7 +299,7 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & integer(I4B), intent(in) :: kiter !< solution outer iteration number ! -- local integer(I4B) :: n - real(DP) :: hhcof, rrhs + real(DP) :: rrhs real(DP) :: vcell real(DP) :: decay_rate ! @@ -310,7 +310,6 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & if (this%ibound(n) <= 0) cycle ! ! -- set variables - hhcof = DZERO rrhs = DZERO vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) ! From f5c2c9404fdc47567e3beb98c539e6422e1f0bd4 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 16:08:57 -0800 Subject: [PATCH 12/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930609444 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 93b3b0a820f..067b93ffe5e 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -67,7 +67,7 @@ module GweEstModule procedure :: est_ar procedure :: est_fc procedure :: est_fc_sto - procedure :: est_fc_dcy + procedure :: est_fc_dcy_water procedure :: est_fc_dcy_solid procedure :: est_cq procedure :: est_cq_sto @@ -178,8 +178,8 @@ subroutine est_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, & ! ! -- decay contribution if (this%idcy /= 0) then - call this%est_fc_dcy(nodes, cold, cnew, nja, matrix_sln, idxglo, & - rhs, kiter) + call this%est_fc_dcy_water(nodes, cold, cnew, nja, matrix_sln, idxglo, & + rhs, kiter) call this%est_fc_dcy_solid(nodes, cold, nja, matrix_sln, idxglo, rhs, & cnew, kiter) end if @@ -237,8 +237,8 @@ end subroutine est_fc_sto !! !! Method to calculate and fill decay coefficients for the package. !< - subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & - idxglo, rhs, kiter) + subroutine est_fc_dcy_water(this, nodes, cold, cnew, nja, matrix_sln, & + idxglo, rhs, kiter) ! -- dummy class(GweEstType) :: this !< GweEstType object integer, intent(in) :: nodes !< number of nodes @@ -279,7 +279,7 @@ subroutine est_fc_dcy(this, nodes, cold, cnew, nja, matrix_sln, & end if ! end do - end subroutine est_fc_dcy + end subroutine est_fc_dcy_water !> @ brief Fill solid decay coefficient method for package !! From 2dc749c20e9c4eae6ccaab045e4136909bae8543 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 16:25:44 -0800 Subject: [PATCH 13/18] Fixed for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930649875 --- doc/mf6io/mf6ivar/dfn/gwe-est.dfn | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn index 6172f361c66..4115e7ac40a 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn @@ -14,7 +14,7 @@ type keyword reader urword optional true longname activate zero-order decay in aqueous phase -description is a text keyword to indicate that zero-order decay will occur in the aqueous phase. Use of this keyword requires that DECAY\_WATER is specified in the GRIDDATA block. +description is a text keyword to indicate that zero-order decay will occur in the aqueous phase. That is, decay occurs in the water and is a rate per volume of water only, not per volume of aquifer (i.e., grid cell). Use of this keyword requires that DECAY\_WATER is specified in the GRIDDATA block. block options name zero_order_decay_solid @@ -22,7 +22,7 @@ type keyword reader urword optional true longname activate zero-order decay in solid phase -description is a text keyword to indicate that zero-order decay will occur in the solid phase. Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block. +description is a text keyword to indicate that zero-order decay will occur in the solid phase. That is, decay occurs in the solid and is a rate per mass of solid only, not per volume of aquifer (i.e., grid cell). Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block. block options name density_water @@ -73,7 +73,7 @@ reader readarray layered true optional true longname aqueous phase decay rate coefficient -description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the aqueous phase is energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block. +description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of zero-order decay in the aqueous phase are energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block. block griddata name decay_solid @@ -83,7 +83,7 @@ reader readarray layered true optional true longname solid phase decay rate coefficient -description is the rate coefficient for zero-order decay for the solid phase. A negative value indicates heat (energy) production. The dimensions of decay for zero-order decay in the solid phase is energy per length cubed per time. Zero-order decay in the solid phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_SOLID is specified in the options block. +description is the rate coefficient for zero-order decay for the solid phase. A negative value indicates heat (energy) production. The dimensions of zero-order decay in the solid phase are energy per mass of solid per time. Zero-order decay in the solid phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_SOLID is specified in the options block. block griddata name heat_capacity_solid From 179b7de19892391555506071ae122dbf26e087bd Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Wed, 29 Jan 2025 21:33:55 -0800 Subject: [PATCH 14/18] New test in response to https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1930582264 --- autotest/test_gwe_decay02.py | 162 +++++++++++++++++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 autotest/test_gwe_decay02.py diff --git a/autotest/test_gwe_decay02.py b/autotest/test_gwe_decay02.py new file mode 100644 index 00000000000..e2c8b4717b3 --- /dev/null +++ b/autotest/test_gwe_decay02.py @@ -0,0 +1,162 @@ +""" +Test problem for decay of energy in EST package of GWE. Compares energy loss +to that of an ESL boundary +""" + +# Imports + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +# Base simulation and model name and workspace + +cases = ["decay"] + +# Model units +length_units = "meters" +time_units = "seconds" + +rho_w = 1000 # kg/m^3 +rho_s = 2500 # kg/m^3 +n = 0.2 # - +gamma_w = -1000 # J/s/m^3, arbitrary value for zero-order aqueous heat production +gamma_s = -0.1 # J/s/kg, arbitrary value for zero-order solid heat production +c_w = 4000 # J/kg/degC +c_s = 1000 # J/kg/decC +T0 = 0 # degC + +nrow = 1 +ncol = 1 +nlay = 1 +delr = 1 # m +delc = 1 # m +top = 1 # m +botm = 0 # m + +perlen = 86400 # s +nstp = 20 + + +def add_gwe(sim, gwename, add_esl=False): + gwe = flopy.mf6.ModflowGwe( + sim, + modelname=gwename, + save_flows=True, + model_nam_file=f"{gwename}.nam", + ) + dis = flopy.mf6.ModflowGwedis( + gwe, nrow=nrow, ncol=ncol, nlay=nlay, delr=delr, delc=delc, top=top, botm=botm + ) + ic = flopy.mf6.ModflowGweic(gwe, strt=T0) + if not add_esl: + zero_order_decay_water = True + zero_order_decay_solid = True + else: + zero_order_decay_water = False + zero_order_decay_solid = False + + esl_spd = { + 0: [[(0, 0, 0), 400]], + } + esl = flopy.mf6.ModflowGweesl( + gwe, + stress_period_data=esl_spd, + pname="ESL", + filename=f"{gwename}.esl", + ) + + est = flopy.mf6.ModflowGweest( + gwe, + zero_order_decay_water=zero_order_decay_water, + zero_order_decay_solid=zero_order_decay_solid, + density_water=rho_w, + density_solid=rho_s, + heat_capacity_water=c_w, + heat_capacity_solid=c_s, + porosity=n, + decay_water=gamma_w, + decay_solid=gamma_s, + ) + + oc = flopy.mf6.ModflowGweoc( + gwe, + budget_filerecord=f"{gwe.name}.bud", + temperature_filerecord=f"{gwe.name}.ucn", + printrecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")], + saverecord=[("BUDGET", "ALL"), ("TEMPERATURE", "ALL")], + ) + + return sim + + +def build_models(idx, test): + # Base MF6 GWF model type + ws = test.workspace + name = cases[idx] + gwename = "gwe-" + name + + print(f"Building MF6 model...{name}") + + sim = flopy.mf6.MFSimulation( + sim_name="heat", + sim_ws=ws, + exe_name="mf6", + version="mf6", + ) + tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(perlen, nstp, 1.0)]) + ims = flopy.mf6.ModflowIms( + sim, complexity="SIMPLE", inner_dvclose=0.001 + ) # T can not become negative in this model + + # add first GWE model + sim = add_gwe(sim, gwename + "-1", add_esl=False) + # add second GWE model + sim = add_gwe(sim, gwename + "-2", add_esl=True) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + ws = test.workspace + + msg = ( + "Differences detected between the simulated results for zeroth-order " + "energy decay and the ESL Package. The respective approaches are " + "expected to give the same answer." + ) + + # read transport results from GWE model + name = cases[idx] + gwename1 = "gwe-" + name + "-1" + gwename2 = "gwe-" + name + "-2" + + # Get the MF6 temperature output + sim = test.sims[0] + gwe1 = sim.get_model(gwename1) + temp_ts1 = gwe1.output.temperature().get_ts((0, 0, 0)) + t1 = temp_ts1[:, 0] + + gwe2 = sim.get_model(gwename2) + temp_ts2 = gwe2.output.temperature().get_ts((0, 0, 0)) + t2 = temp_ts2[:, 0] + + assert np.isclose(temp_ts1[-1, 1], temp_ts2[-1, 1], atol=1e-10), msg + + +# - No need to change any code below +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() From 495132eb0af4306ef677b2e0e9a98119b0b5b034 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Thu, 30 Jan 2025 08:47:49 -0800 Subject: [PATCH 15/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1935834832 --- autotest/test_gwe_decay02.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/autotest/test_gwe_decay02.py b/autotest/test_gwe_decay02.py index e2c8b4717b3..a21ce12b481 100644 --- a/autotest/test_gwe_decay02.py +++ b/autotest/test_gwe_decay02.py @@ -57,8 +57,9 @@ def add_gwe(sim, gwename, add_esl=False): zero_order_decay_water = False zero_order_decay_solid = False + esl_amt = n * gamma_w + (1 - n) * gamma_s * rho_s esl_spd = { - 0: [[(0, 0, 0), 400]], + 0: [[(0, 0, 0), abs(esl_amt)]], } esl = flopy.mf6.ModflowGweesl( gwe, From c52ed0ea71b13166e9c5ad45a5c077311d2822fa Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Thu, 30 Jan 2025 09:08:33 -0800 Subject: [PATCH 16/18] Adjustments applied in response to https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1935842327 and https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1935846602 --- doc/mf6io/mf6ivar/dfn/gwe-est.dfn | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn index 4115e7ac40a..c975ab0d424 100644 --- a/doc/mf6io/mf6ivar/dfn/gwe-est.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwe-est.dfn @@ -22,7 +22,7 @@ type keyword reader urword optional true longname activate zero-order decay in solid phase -description is a text keyword to indicate that zero-order decay will occur in the solid phase. That is, decay occurs in the solid and is a rate per mass of solid only, not per volume of aquifer (i.e., grid cell). Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block. +description is a text keyword to indicate that zero-order decay will occur in the solid phase. That is, decay occurs in the solid and is a rate per mass (not volume) of solid only. Use of this keyword requires that DECAY\_SOLID is specified in the GRIDDATA block. block options name density_water @@ -73,7 +73,7 @@ reader readarray layered true optional true longname aqueous phase decay rate coefficient -description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of zero-order decay in the aqueous phase are energy per length cubed per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block. +description is the rate coefficient for zero-order decay for the aqueous phase of the mobile domain. A negative value indicates heat (energy) production. The dimensions of zero-order decay in the aqueous phase are energy per length cubed (volume of water) per time. Zero-order decay in the aqueous phase will have no effect on simulation results unless ZERO\_ORDER\_DECAY\_WATER is specified in the options block. block griddata name decay_solid From eb62e823340d4f5cccf0ad15f2a0fd20c9a55dd6 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Thu, 30 Jan 2025 09:22:40 -0800 Subject: [PATCH 17/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1935989639 --- autotest/test_gwe_decay02.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/autotest/test_gwe_decay02.py b/autotest/test_gwe_decay02.py index a21ce12b481..35e2e2c12dd 100644 --- a/autotest/test_gwe_decay02.py +++ b/autotest/test_gwe_decay02.py @@ -59,7 +59,7 @@ def add_gwe(sim, gwename, add_esl=False): esl_amt = n * gamma_w + (1 - n) * gamma_s * rho_s esl_spd = { - 0: [[(0, 0, 0), abs(esl_amt)]], + 0: [[(0, 0, 0), -esl_amt]], } esl = flopy.mf6.ModflowGweesl( gwe, From b7a233069dbd2e2a7c809059d8e558e4ac5d2d21 Mon Sep 17 00:00:00 2001 From: emorway-usgs Date: Fri, 31 Jan 2025 08:51:38 -0800 Subject: [PATCH 18/18] Fix for review comment https://github.com/MODFLOW-USGS/modflow6/pull/2155#discussion_r1937356433 --- src/Model/GroundWaterEnergy/gwe-est.f90 | 73 ++++++++++++++----------- 1 file changed, 42 insertions(+), 31 deletions(-) diff --git a/src/Model/GroundWaterEnergy/gwe-est.f90 b/src/Model/GroundWaterEnergy/gwe-est.f90 index 067b93ffe5e..cd5e4356e62 100644 --- a/src/Model/GroundWaterEnergy/gwe-est.f90 +++ b/src/Model/GroundWaterEnergy/gwe-est.f90 @@ -31,6 +31,16 @@ module GweEstModule character(len=LENBUDTXT), dimension(NBDITEMS) :: budtxt data budtxt/' STORAGE-CELLBLK', ' DECAY-AQUEOUS', ' DECAY-SOLID'/ + !> @brief Enumerator that defines the decay options + !< + ENUM, BIND(C) + ENUMERATOR :: DECAY_OFF = 0 !< Decay (or production) of thermal energy inactive (default) + ENUMERATOR :: DECAY_ZERO_ORDER = 2 !< Zeroth-order decay + ENUMERATOR :: DECAY_WATER = 1 !< Zeroth-order decay in water only + ENUMERATOR :: DECAY_SOLID = 2 !< Zeroth-order decay in solid only + ENUMERATOR :: DECAY_BOTH = 3 !< Zeroth-order decay in water and solid + END ENUM + !> @ brief Energy storage and transfer !! !! Data and methods for handling changes in temperature @@ -177,7 +187,7 @@ subroutine est_fc(this, nodes, cold, nja, matrix_sln, idxglo, cnew, & call this%est_fc_sto(nodes, cold, nja, matrix_sln, idxglo, rhs) ! ! -- decay contribution - if (this%idcy /= 0) then + if (this%idcy == DECAY_ZERO_ORDER) then call this%est_fc_dcy_water(nodes, cold, cnew, nja, matrix_sln, idxglo, & rhs, kiter) call this%est_fc_dcy_solid(nodes, cold, nja, matrix_sln, idxglo, rhs, & @@ -267,8 +277,8 @@ subroutine est_fc_dcy_water(this, nodes, cold, cnew, nja, matrix_sln, & swtpdt = this%fmi%gwfsat(n) ! ! -- add zero-order decay rate terms to accumulators - if (this%idcy == 2 .and. (this%idcysrc == 1 .or. & - this%idcysrc == 3)) then + if (this%idcy == DECAY_ZERO_ORDER .and. (this%idcysrc == DECAY_WATER .or. & + this%idcysrc == DECAY_BOTH)) then ! decay_rate = this%decay_water(n) ! -- This term does get divided by eqnsclfac for fc purposes because it @@ -314,8 +324,8 @@ subroutine est_fc_dcy_solid(this, nodes, cold, nja, matrix_sln, idxglo, & vcell = this%dis%area(n) * (this%dis%top(n) - this%dis%bot(n)) ! ! -- account for zero-order decay rate terms in rhs - if (this%idcy == 2 .and. (this%idcysrc == 2 .or. & - this%idcysrc == 3)) then + if (this%idcy == DECAY_ZERO_ORDER .and. (this%idcysrc == DECAY_SOLID .or. & + this%idcysrc == DECAY_BOTH)) then ! ! -- negative temps are currently not checked for or prevented since a ! user can define a temperature scale of their own choosing. if @@ -347,11 +357,11 @@ subroutine est_cq(this, nodes, cnew, cold, flowja) call this%est_cq_sto(nodes, cnew, cold, flowja) ! ! -- decay - if (this%idcy /= 0) then - if (this%idcysrc == 1 .or. this%idcysrc == 3) then + if (this%idcy == DECAY_ZERO_ORDER) then + if (this%idcysrc == DECAY_WATER .or. this%idcysrc == DECAY_BOTH) then call this%est_cq_dcy(nodes, cnew, cold, flowja) end if - if (this%idcysrc == 2 .or. this%idcysrc == 3) then + if (this%idcysrc == DECAY_SOLID .or. this%idcysrc == DECAY_BOTH) then call this%est_cq_dcy_solid(nodes, cnew, cold, flowja) end if end if @@ -445,8 +455,8 @@ subroutine est_cq_dcy(this, nodes, cnew, cold, flowja) hhcof = DZERO rrhs = DZERO ! -- zero order decay aqueous phase - if (this%idcy == 2 .and. & - (this%idcysrc == 1 .or. this%idcysrc == 3)) then + if (this%idcy == DECAY_ZERO_ORDER .and. & + (this%idcysrc == DECAY_WATER .or. this%idcysrc == DECAY_BOTH)) then decay_rate = this%decay_water(n) ! -- this term does NOT get multiplied by eqnsclfac for cq purposes ! because it should already be a rate of energy @@ -494,8 +504,8 @@ subroutine est_cq_dcy_solid(this, nodes, cnew, cold, flowja) hhcof = DZERO rrhs = DZERO ! -- first-order decay (idcy=1) is not supported for temperature modeling - if (this%idcy == 2 .and. & - (this%idcysrc == 2 .or. this%idcysrc == 3)) then ! zero order decay in the solid phase + if (this%idcy == DECAY_ZERO_ORDER .and. & + (this%idcysrc == DECAY_SOLID .or. this%idcysrc == DECAY_BOTH)) then ! zero order decay in the solid phase decay_rate = this%decay_solid(n) ! -- this term does NOT get multiplied by eqnsclfac for cq purposes ! because it should already be a rate of energy @@ -531,14 +541,14 @@ subroutine est_bd(this, isuppress_output, model_budget) isuppress_output, rowlabel=this%packName) ! ! -- dcy - if (this%idcy /= 0) then - if (this%idcysrc == 1 .or. this%idcysrc == 3) then + if (this%idcy == DECAY_ZERO_ORDER) then + if (this%idcysrc == DECAY_WATER .or. this%idcysrc == DECAY_BOTH) then ! -- aqueous phase call rate_accumulator(this%ratedcyw, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(2), & isuppress_output, rowlabel=this%packName) end if - if (this%idcysrc == 2 .or. this%idcysrc == 3) then + if (this%idcysrc == DECAY_SOLID .or. this%idcysrc == DECAY_BOTH) then ! -- solid phase call rate_accumulator(this%ratedcys, rin, rout) call model_budget%addentry(rin, rout, delt, budtxt(3), & @@ -583,14 +593,14 @@ subroutine est_ot_flow(this, icbcfl, icbcun) nwidthp, editdesc, dinact) ! ! -- dcy - if (this%idcy /= 0) then - if (this%idcysrc == 1 .or. this%idcysrc == 3) then + if (this%idcy == DECAY_ZERO_ORDER) then + if (this%idcysrc == DECAY_WATER .or. this%idcysrc == DECAY_BOTH) then ! -- aqueous phase call this%dis%record_array(this%ratedcyw, this%iout, iprint, & -ibinun, budtxt(2), cdatafmp, nvaluesp, & nwidthp, editdesc, dinact) end if - if (this%idcysrc == 2 .or. this%idcysrc == 3) then + if (this%idcysrc == DECAY_SOLID .or. this%idcysrc == DECAY_BOTH) then ! -- solid phase call this%dis%record_array(this%ratedcys, this%iout, iprint, & -ibinun, budtxt(3), cdatafmp, nvaluesp, & @@ -687,7 +697,7 @@ subroutine allocate_arrays(this, nodes) call mem_allocate(this%rhos, nodes, 'RHOS', this%memoryPath) ! ! -- dcy - if (this%idcy == 0) then + if (this%idcy == DECAY_OFF) then call mem_allocate(this%ratedcyw, 1, 'RATEDCYW', this%memoryPath) call mem_allocate(this%ratedcys, 1, 'RATEDCYS', this%memoryPath) call mem_allocate(this%decay_water, 1, 'DECAY_WATER', this%memoryPath) @@ -763,23 +773,24 @@ subroutine read_options(this) case ('SAVE_FLOWS') this%ipakcb = -1 write (this%iout, fmtisvflow) - case ('FIRST_ORDER_DECAY') - this%idcy = 1 - write (this%iout, fmtidcy1) case ('ZERO_ORDER_DECAY_WATER') - this%idcy = 2 + this%idcy = DECAY_ZERO_ORDER + ! -- idcysrc > 0 indicates decay in the solid phase is active + ! in which case the idcysrc should now be upgraded to both phases if (this%idcysrc > IZERO) then - this%idcysrc = 3 + this%idcysrc = DECAY_BOTH else - this%idcysrc = 1 + this%idcysrc = DECAY_WATER end if write (this%iout, fmtidcy2) case ('ZERO_ORDER_DECAY_SOLID') - this%idcy = 2 + this%idcy = DECAY_ZERO_ORDER + ! -- idcysrc > 0 indicates decay in active in water in which case + ! the idcysrc should now be upgraded to both phases if (this%idcysrc > IZERO) then - this%idcysrc = 3 + this%idcysrc = DECAY_BOTH else - this%idcysrc = 2 + this%idcysrc = DECAY_SOLID end if write (this%iout, fmtidcy3) case ('HEAT_CAPACITY_WATER') @@ -867,7 +878,7 @@ subroutine read_data(this) aname(1)) lname(1) = .true. case ('DECAY_WATER') - if (this%idcy == 0) & + if (this%idcy == DECAY_OFF) & call mem_reallocate(this%decay_water, this%dis%nodes, 'DECAY_WATER', & trim(this%memoryPath)) call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & @@ -875,7 +886,7 @@ subroutine read_data(this) aname(2)) lname(2) = .true. case ('DECAY_SOLID') - if (this%idcy == 0) & + if (this%idcy == DECAY_OFF) & call mem_reallocate(this%decay_solid, this%dis%nodes, 'DECAY_SOLID', & trim(this%memoryPath)) call this%dis%read_grid_array(line, lloc, istart, istop, this%iout, & @@ -920,7 +931,7 @@ subroutine read_data(this) end if ! ! -- Check for required decay/production rate coefficients - if (this%idcy > 0) then + if (this%idcy == DECAY_ZERO_ORDER) then if (.not. lname(2) .and. .not. lname(3)) then write (errmsg, '(a)') 'Zero order decay in either the aqueous & &or solid phase is active but the corresponding zero-order &