From 3a6b1ef7988184f5e00ae7ae9c00640c92cfad06 Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Thu, 8 Feb 2024 11:49:51 -0800 Subject: [PATCH] Adding multi-aquifer well energy transport (MWE) package. Includes new autotest --- autotest/test_gwe_mwe_conduction.py | 570 +++++++++ doc/Common/gwe-mweobs.tex | 16 +- doc/mf6io/gwe/gwe.tex | 4 + doc/mf6io/gwe/mwe.tex | 55 + doc/mf6io/gwe/namefile.tex | 1 + doc/mf6io/mf6ivar/dfn/gwe-mwe.dfn | 447 ++++++++ .../mf6ivar/examples/gwe-mwe-example-obs.dat | 43 + .../mf6ivar/examples/gwe-mwe-example.dat | 24 + doc/mf6io/mf6ivar/md/mf6ivar.md | 35 + doc/mf6io/mf6ivar/mf6ivar.py | 1 + doc/mf6io/mf6ivar/tex/appendixA.tex | 4 + doc/mf6io/mf6ivar/tex/gwe-mwe-desc.tex | 92 ++ doc/mf6io/mf6ivar/tex/gwe-mwe-options.dat | 15 + doc/mf6io/mf6ivar/tex/gwe-mwe-packagedata.dat | 5 + doc/mf6io/mf6ivar/tex/gwe-mwe-period.dat | 5 + doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex | 2 +- make/makefile | 1 + msvs/mf6core.vfproj | 1 + src/Model/GroundWaterEnergy/gwe1.f90 | 10 +- src/Model/GroundWaterEnergy/gwe1mwe1.f90 | 1018 +++++++++++++++++ src/Model/TransportModel/tsp1apt1.f90 | 4 +- src/meson.build | 1 + 22 files changed, 2338 insertions(+), 16 deletions(-) create mode 100644 autotest/test_gwe_mwe_conduction.py create mode 100644 doc/mf6io/gwe/mwe.tex create mode 100644 doc/mf6io/mf6ivar/dfn/gwe-mwe.dfn create mode 100644 doc/mf6io/mf6ivar/examples/gwe-mwe-example-obs.dat create mode 100644 doc/mf6io/mf6ivar/examples/gwe-mwe-example.dat create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mwe-desc.tex create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mwe-options.dat create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mwe-packagedata.dat create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mwe-period.dat create mode 100644 src/Model/GroundWaterEnergy/gwe1mwe1.f90 diff --git a/autotest/test_gwe_mwe_conduction.py b/autotest/test_gwe_mwe_conduction.py new file mode 100644 index 00000000000..b5570a0527f --- /dev/null +++ b/autotest/test_gwe_mwe_conduction.py @@ -0,0 +1,570 @@ +# Test mwe package. Looks at wetted area of well for calculating +# heat conduction exchange. This test is related to test_gwf_maw06.py but +# with some further customization for testing GWE capabilities. +# - Test has 0 flow conductance with gw cells in layers 1-3 +# - Test uses MAW to inject water into layer 4 (bottom layer) +# - Test includes conductive exchng only between MWE feature and +# layers 1 to 3 +# - Water extracted by normal WEL features in bottom corners of model +# - Water table saturates only ~1/2 of top layer; therefore +# conductive exchange + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework + + +def process_line(line): + m_arr = line.strip().split() + if any("=" in itm and len(itm) > 1 for itm in m_arr): + m_arr = [ + float(itm.split("=")[-1]) if len(itm.split("=")) > 1 else itm + for itm in m_arr + ] + nm = m_arr[-2] + else: + nm = m_arr[-3] + val = m_arr[-1] + return {nm: float(val)} + + +def get_bud(fname, srchStr): + in_bud_lst = {} + out_bud_lst = {} + with open(fname, "r") as f: + for line in f: + if srchStr in line: + # Read the package budget + line = next(f) + while not "TOTAL IN =" in line: + if "=" in line: + in_bud_lst.update(process_line(line)) + + line = next(f) + + # Get "total in" + dct = process_line(line) + T_in = dct["IN"] + + line = next(f) + while not "TOTAL OUT =" in line: + if "=" in line: + out_bud_lst.update(process_line(line)) + + line = next(f) + + # Get "total out" + dct = process_line(line) + T_out = dct["OUT"] + + break + + return T_in, T_out, in_bud_lst, out_bud_lst + + +def get_welbore_heat_flow(fname, srchStr): + ener_Q = [] + with open(fname, "r") as f: + for line in f: + if srchStr in line: + # Read an established format + for i in np.arange(3): # read & discard 3 lines + line = next(f) + for i in np.arange( + 4 + ): # read & digest 4 lines of needed output + line = next(f) + m_arr = line.strip().split() + ener_Q.append([int(m_arr[0]), float(m_arr[2])]) + break + + return np.array(ener_Q) + + +def trenddetector(list_of_index, array_of_data, order=1): + result = np.polyfit(list_of_index, list(array_of_data), order) + slope = result[-2] + return float(slope) + + +cases = ["mwe_01"] +mawstrt = 3.5 + +# Flow related parameters +lx = 70.0 +ly = 70.0 +nlay = 4 +nrow = 7 +ncol = 7 +nper = 1 +delc = ly / nrow +delr = lx / ncol +top = 4.0 +botm = [3.0, 2.0, 1.0, 0.0] +strt = 3.5 +transient = {0: True} + +perlen = [10.0] +nstp = [10] +tsmult = [1.0] + +Kh = [1.0, 1.0, 1e-6, 100] +Kv = [1.0, 1.0, 1e-6, 100] + +Sy = 0.3 +Ss = 0.0 + +# Transport related parameters +mixelm = 0 # Upstream vs TVD (Upstream selected) +strttemp = 1.0 # Initial temperature ($^{\circ}C$) +porosity = Sy # porosity (unitless) +ktw = 0.5918 # Thermal Conductivity of Water ($W/m/^{\circ}C$) +kts = 0.2700 # Thermal Conductivity of Aquifer Solids ($W/m/^{\circ}C$) +rhow = 1000 # Density of water ($kg/m^3$) +rhos = 2650 # Density of the aquifer material ($kg/m^3$) +Cpw = 4180 # Heat Capacity of water ($J/kg/C$) +Cps = 880 # Heat capacity of the solids ($J/kg/C$) +lhv = 2454000.0 # Latent heat of vaporization ($J/kg$) +K_therm_maw = 1.5 # Thermal conductivity of the lakebed material ($W/m/C$) +wthkcnd = 0.01 +mawradius = 0.1 +mawbottom = 0.0 + +# Solver settings +nouter, ninner = 700, 10 +hclose, rclose, relax = 1e-8, 1e-6, 0.97 + + +def build_models(idx, test): + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = cases[idx] + + # Instantiate MODFLOW 6 simulation + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # Instantiate Time Discretization package + flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # Instantiate Groundwater Flow model + gwfname = "gwf-" + name + gwename = "gwe-" + name + newtonoptions = "NEWTON UNDER_RELAXATION" + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, newtonoptions=newtonoptions + ) + + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="SIMPLE", + under_relaxation_gamma=0.1, + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwfname}.ims", + ) + sim.register_ims_package(ims, [gwfname]) + + # Instantiate Discretization package + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + length_units="METERS", + pname="DIS", + filename=f"{gwfname}.dis", + ) + + # Instantiate Initial Conditions package + flopy.mf6.ModflowGwfic( + gwf, + strt=strt, + pname="IC", + filename=f"{gwfname}.ic", + ) + + # Instantiate Node Property Flow package + flopy.mf6.ModflowGwfnpf( + gwf, + xt3doptions=False, + save_flows=True, + save_specific_discharge=True, + icelltype=1, + k=Kh, + k33=Kv, + pname="NPF", + filename=f"{gwfname}.npf", + ) + + # Instantiate Storage package + flopy.mf6.ModflowGwfsto( + gwf, + sy=Sy, + ss=Ss, + iconvert=1, + transient=transient, + pname="STO", + filename=f"{gwfname}.sto", + ) + + # Instantiate Well package (for extracting MAW-injected water) + wellist = [] + for i in [0, nrow - 1]: + for j in [0, ncol - 1]: + wellist.append([(nlay - 1, i, j), -7.0, 0.01]) + + flopy.mf6.ModflowGwfwel( + gwf, + stress_period_data=wellist, + print_input=True, + print_flows=True, + save_flows=False, + pname="WEL", + auxiliary="TEMPERATURE", + filename=f"{gwfname}.wel", + ) + + # Instantiate Multi-Aquifer Well package (injects water) + mstrt = mawstrt + mawcondeqn = "SPECIFIED" + mawngwfnodes = nlay + # + mawpackagedata = [ + [0, mawradius, mawbottom, mstrt, mawcondeqn, mawngwfnodes] + ] + # + conncond = [0.0, 0.0, 0.0, 1000.0] + mawconnectiondata = [ + [0, icon, (icon, 3, 3), top, mawbottom, conncond[icon], -999.0] + for icon in range(nlay) + ] + # + mawperioddata = [[0, "STATUS", "ACTIVE"], [0, "RATE", 28.0]] + maw = flopy.mf6.ModflowGwfmaw( + gwf, + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + head_filerecord=f"{gwfname}.maw.bin", + budget_filerecord=f"{gwfname}.maw.cbc", + packagedata=mawpackagedata, + connectiondata=mawconnectiondata, + perioddata=mawperioddata, + pname="MAW-1", + ) + opth = f"{gwfname}.maw.obs" + obsdata = { + f"{gwfname}.maw.obs.csv": [ + ("whead", "head", (0,)), + ] + } + maw.obs.initialize( + filename=opth, digits=20, print_input=True, continuous=obsdata + ) + + # Instantiate Output Control package + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{gwfname}.cbc", + head_filerecord=f"{gwfname}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[ + ( + "HEAD", + "ALL", + ), + ( + "BUDGET", + "ALL", + ), + ], + printrecord=[ + ( + "HEAD", + "ALL", + ), + ( + "BUDGET", + "ALL", + ), + ], + ) + + # Create GWE model + # ---------------- + gwe = flopy.mf6.ModflowGwe( + sim, modelname=gwename, model_nam_file="{}.nam".format(gwename) + ) + gwe.name_file.save_flows = True + + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + filename=f"{gwename}.ims", + ) + sim.register_ims_package(imsgwe, [gwename]) + + # Instantiating MODFLOW 6 enregy transport discretization package + flopy.mf6.ModflowGwedis( + gwe, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename="{}.dis".format(gwename), + ) + + # Instantiating MODFLOW 6 energy transport initial temperatures + flopy.mf6.ModflowGweic( + gwe, strt=strttemp, filename="{}.ic".format(gwename) + ) + + # Instantiate mobile storage and transfer package + flopy.mf6.ModflowGweest( + gwe, + porosity=porosity, + cps=Cps, + rhos=rhos, + packagedata=[Cpw, rhow, lhv], + pname="MST-1", + filename=f"{gwename}.mst", + ) + + # Instantiating MODFLOW 6 energy transport advection package + if mixelm == 0: + scheme = "UPSTREAM" + elif mixelm == -1: + scheme = "TVD" + else: + raise Exception() + + # Instantiate advection package + flopy.mf6.ModflowGweadv( + gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename) + ) + + # Instantiate dispersion package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=True, + ktw=0.5918, + kts=0.2700, + filename="{}.dsp".format(gwename), + ) + + # Instantiate source/sink mixing package + sourcerecarray = [ + ("WEL", "AUX", "TEMPERATURE"), + ] + flopy.mf6.ModflowGwessm( + gwe, sources=sourcerecarray, filename=f"{gwename}.ssm" + ) + + # Instantiating MODFLOW 6 transport output control package + flopy.mf6.ModflowGweoc( + gwe, + budget_filerecord="{}.cbc".format(gwename), + temperature_filerecord="{}.ucn".format(gwename), + temperatureprintrecord=[ + ("COLUMNS", 17, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename="{}.oc".format(gwename), + ) + + # Instantiating MODFLOW 6 multi-well energy transport (mwe) package + # ,, , , + mwepackagedata = [(0, 1.0, K_therm_maw, wthkcnd, "well1")] + + mweperioddata = {0: [(0, "RATE", 40.0)]} + + # note: for specifying lake number, use fortran indexing! + mwe_obs = { + "{}.mweobs".format(gwename): [ + ("MweTemp", "temperature", 1), + ] + } + + flopy.mf6.ModflowGwemwe( + gwe, + flow_package_name="MAW-1", + budget_filerecord=gwename + ".mwe.bud", + boundnames=True, + save_flows=True, + print_input=True, + print_flows=True, + print_temperature=True, + packagedata=mwepackagedata, + mweperioddata=mweperioddata, + observations=mwe_obs, + pname="MWE-1", + filename="{}.mwe".format(gwename), + ) + + # Instantiate GWF-GWE exchange + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename=f"{name}.gwfgwe", + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating results...") + + top_local = [4.0, 3.0, 2.0, 1.0] + botm_local = [3.0, 2.0, 1.0, 0.0] + + # calculate volume of water and make sure it is conserved + name = cases[idx] + gwfname = "gwf-" + name + fname = gwfname + ".maw.bin" + fname = os.path.join(test.workspace, fname) + assert os.path.isfile(fname) + bobj = flopy.utils.HeadFile(fname, text="HEAD") + stage = bobj.get_alldata().flatten() + + name = cases[idx] + gwfname = "gwf-" + name + fname = gwfname + ".maw.cbc" + fname = os.path.join(test.workspace, fname) + assert os.path.isfile(fname) + bobj = flopy.utils.CellBudgetFile(fname, precision="double") + gwfarea = bobj.get_data(text="GWF") + + # Retrieve simulated temperature for the multi-aquifer well + gwename = "gwe-" + name + fname = gwename + ".mweobs" + mwtemp_file = os.path.join(test.workspace, fname) + assert os.path.isfile(mwtemp_file) + mwtemp = np.genfromtxt(mwtemp_file, names=True, delimiter=",") + mwtemp = mwtemp["MWETEMP"].astype(float).reshape((mwtemp.size, 1)) + + # Retrieve gw temperatures + fname = gwename + ".ucn" + fname = os.path.join(test.workspace, fname) + assert os.path.isfile(fname) + gwtempobj = flopy.utils.HeadFile( + fname, precision="double", text="TEMPERATURE" + ) + gwe_temps = gwtempobj.get_alldata() + + # Calculate conductive exchange external to MF6 and compare to MF6 values + # Evaluates first time step only + wellbore_cnd_time1 = [] + for i in np.arange(nlay): + if stage[0] > top_local[i]: + thk = top_local[i] - botm_local[i] + elif stage[0] > botm_local[i]: + thk = stage[0] - botm_local[i] + else: + thk = 0 + + # Check that MF6 (GWF) wellbore wetted area matches explicitly calc + + wa = 2 * mawradius * np.pi * thk + welborecnd = K_therm_maw * wa / wthkcnd + gw_temp = gwe_temps[0, i, 3, 3] + deltaT = mwtemp[0][0] - gw_temp + + wellbore_cnd_time1.append(welborecnd * deltaT) + + # Retrieve budget + fname = os.path.join(test.workspace, gwename + ".lst") + srchStr = ( + "MWE-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 1, " + "STRESS PERIOD 1" + ) + T_in, T_out, in_bud_lst, out_bud_lst = get_bud(fname, srchStr) + assert np.isclose( + T_in, T_out, atol=0.1 + ), "There is a heat budget discrepancy where there shouldn't be" + + msg1 = ( + "Conductive heat exchanges calculated explicitly and by MF6 " + "do not match" + ) + msg2 = ( + "Individually summing well bore 'heat flows' is not matching " + "the global budget heat flow into the aquifer" + ) + msg3 = "Groundwater should be warming, but isn't" + + # Get MF6 saved wellbore heat "flows" + srchStr = "MWE PACKAGE (MWE-1) FLOW RATES PERIOD 1 STEP 1" + wbcnd_mf6 = get_welbore_heat_flow(fname, srchStr) + + # Check top 3 layers (4th layer handled different) + for i in np.arange(nlay - 1): + assert np.isclose( + wbcnd_mf6[i, 1], round(wellbore_cnd_time1[i], 4) + ), msg1 + + # Layer 4 "heat flow" includes convection and conduction, compare + # "heat flow" from all layers to global budget line item 'IN: GWF' + glob_bud_gw_in = out_bud_lst["GWF"] + out_bud_lst["WELLBORE-COND"] + mwe_out = wbcnd_mf6.sum(axis=0)[1] + assert np.isclose(mwe_out, glob_bud_gw_in, rtol=1e-9), msg2 + + # Ensure that temperatures near the injection point are rising + slp = trenddetector(np.arange(gwe_temps.shape[0]), gwe_temps[:, 3, 3, 3]) + assert slp > 0.0, msg3 + + +# - 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/Common/gwe-mweobs.tex b/doc/Common/gwe-mweobs.tex index 3594db046cf..a2c51bf6a3b 100644 --- a/doc/Common/gwe-mweobs.tex +++ b/doc/Common/gwe-mweobs.tex @@ -1,13 +1,13 @@ % general APT observations -MWE & temperature & ifno or boundname & -- & Well temperature. If boundname is specified, boundname must be unique for each well. \\ +MWE & temperature & mawno or boundname & -- & Well temperature. If boundname is specified, boundname must be unique for each well. \\ %flowjaface not included -MWE & storage & ifno or boundname & -- & Simulated energy storage flow rate for a well or group of wells. \\ -MWE & constant & ifno or boundname & -- & Simulated energy constant-flow rate for a well or group of wells. \\ -MWE & from-mvr & ifno or boundname & -- & Simulated energy inflow into a well or group of wells from the MVE package. Energy inflow is calculated as the product of provider temperature and the mover flow rate. \\ -MWE & mwt & ifno or boundname & \texttt{iconn} or -- & Energy flow rate for a well or group of wells and its aquifer connection(s). If boundname is not specified for ID, then the simulated well-aquifer flow rate at a specific well connection is observed. In this case, ID2 must be specified and is the connection number \texttt{iconn} for well \texttt{ifno}. \\ +MWE & storage & mawno or boundname & -- & Simulated energy storage flow rate for a well or group of wells. \\ +MWE & constant & mawno or boundname & -- & Simulated energy constant-flow rate for a well or group of wells. \\ +MWE & from-mvr & mawno or boundname & -- & Simulated energy inflow into a well or group of wells from the MVE package. Energy inflow is calculated as the product of provider temperature and the mover flow rate. \\ +MWE & mwe & mawno or boundname & \texttt{iconn} or -- & Energy flow rate for a well or group of wells and its aquifer connection(s). If boundname is not specified for ID, then the simulated well-aquifer flow rate at a specific well connection is observed. In this case, ID2 must be specified and is the connection number \texttt{iconn} for well \texttt{mawno}. \\ -% observations specific to the mwt package -MWE & rate & ifno or boundname & -- & Simulated energy flow rate for a well or group of wells. \\ -MWE & fw-rate & ifno or boundname & -- & Simulated energy flow rate for a flowing well or group of flowing wells. \\ +% observations specific to the mwe package +MWE & rate & mawno or boundname & -- & Simulated energy flow rate for a well or group of wells. \\ +MWE & fw-rate & mawno or boundname & -- & Simulated energy flow rate for a flowing well or group of flowing wells. \\ MWE & rate-to-mvr & well or boundname & -- & Simulated energy flow rate that is sent to the MVE Package for a well or group of wells.\\ MWE & fw-rate-to-mvr & well or boundname & -- & Simulated energy flow rate that is sent to the MVE Package from a flowing well or group of flowing wells. \\ diff --git a/doc/mf6io/gwe/gwe.tex b/doc/mf6io/gwe/gwe.tex index 6840cf3baa1..0e7e30211bc 100644 --- a/doc/mf6io/gwe/gwe.tex +++ b/doc/mf6io/gwe/gwe.tex @@ -133,6 +133,10 @@ \subsection{Streamflow Energy Transport (SFE) Package} \subsection{Lake Energy Transport (LKE) Package} \input{gwe/lke} +\newpage +\subsection{Multi-Aquifer Well Energy Transport (MWE) Package} +\input{gwe/mwe} + \newpage \subsection{Flow Model Interface (FMI) Package} \input{gwe/fmi} diff --git a/doc/mf6io/gwe/mwe.tex b/doc/mf6io/gwe/mwe.tex new file mode 100644 index 00000000000..70fc673987d --- /dev/null +++ b/doc/mf6io/gwe/mwe.tex @@ -0,0 +1,55 @@ +Multi-Aquifer Well Energy Transport (MWE) Package information is read from the file that is specified by ``MWE6'' as the file type. There can be as many MWE Packages as necessary for a GWE model. Each MWE Package is designed to work with flows from a corresponding GWF MAW Package. By default \mf uses the MWE package name to determine which MAW Package corresponds to the MWE Package. Therefore, the package name of the MWE Package (as specified in the GWE name file) must match with the name of the corresponding MAW Package (as specified in the GWF name file). Alternatively, the name of the flow package can be specified using the FLOW\_PACKAGE\_NAME keyword in the options block. The GWE MWE Package cannot be used without a corresponding GWF MAW Package. + +The MWE Package does not have a dimensions block; instead, dimensions for the MWE Package are set using the dimensions from the corresponding MAW Package. For example, the MAW Package requires specification of the number of wells (NMAWWELLS). MWE sets the number of wells equal to NMAWWELLS. Therefore, the PACKAGEDATA block below must have NMAWWELLS entries in it. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-mwe-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-mwe-packagedata.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-mwe-period.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/gwe-mwe-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwe-mwe-example.dat} + +\vspace{5mm} +\subsubsection{Available observation types} +Multi-Aquifer Well Energy Transport Package observations include well temperature and all of the terms that contribute to the continuity equation for each well. Additional MWE Package observations include energy flow rates for individual wells, or groups of wells; the well volume (\texttt{volume}); and the conductance for a well-aquifer connection conductance (\texttt{conductance}). The data required for each MWE Package observation type is defined in table~\ref{table:gwe-mweobstype}. Negative and positive values for \texttt{mwe} observations represent a loss from and gain to the GWE model, respectively. For all other flow terms, negative and positive values represent a loss from and gain from the MWE package, respectively. + +\begin{longtable}{p{2cm} p{2.75cm} p{2cm} p{1.25cm} p{7cm}} +\caption{Available MWE Package observation types} \tabularnewline + +\hline +\hline +\textbf{Stress Package} & \textbf{Observation type} & \textbf{ID} & \textbf{ID2} & \textbf{Description} \\ +\hline +\endfirsthead + +\captionsetup{textformat=simple} +\caption*{\textbf{Table \arabic{table}.}{\quad}Available MWE Package observation types.---Continued} \tabularnewline + +\hline +\hline +\textbf{Stress Package} & \textbf{Observation type} & \textbf{ID} & \textbf{ID2} & \textbf{Description} \\ +\hline +\endhead + + +\hline +\endfoot + +\input{../Common/gwe-mweobs.tex} +\label{table:gwe-mweobstype} +\end{longtable} + +\vspace{5mm} +\subsubsection{Example Observation Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwe-mwe-example-obs.dat} + + diff --git a/doc/mf6io/gwe/namefile.tex b/doc/mf6io/gwe/namefile.tex index 12afbd6b30e..4419af99530 100644 --- a/doc/mf6io/gwe/namefile.tex +++ b/doc/mf6io/gwe/namefile.tex @@ -36,6 +36,7 @@ \subsubsection{Explanation of Variables} ESL6 & Energy Source Loading Package & * \\ SFE6 & Streamflow Energy Transport Package & * \\ LKE6 & Lake Energy Transport Package & * \\ +MWE6 & Multi-Aquifer Well Energy Transport Package & * \\ OBS6 & Observations Option \\ \hline \end{tabular*} diff --git a/doc/mf6io/mf6ivar/dfn/gwe-mwe.dfn b/doc/mf6io/mf6ivar/dfn/gwe-mwe.dfn new file mode 100644 index 00000000000..c805b6533fe --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/gwe-mwe.dfn @@ -0,0 +1,447 @@ +# --------------------- gwe mwe options --------------------- +# flopy multi-package + +block options +name flow_package_name +type string +shape +reader urword +optional true +longname keyword to specify name of corresponding flow package +description keyword to specify the name of the corresponding flow package. If not specified, then the corresponding flow package must have the same name as this advanced transport package (the name associated with this package in the GWE name file). + +block options +name auxiliary +type string +shape (naux) +reader urword +optional true +longname keyword to specify aux variables +description REPLACE auxnames {'{#1}': 'Groundwater Energy Transport'} + +block options +name flow_package_auxiliary_name +type string +shape +reader urword +optional true +longname keyword to specify name of temperature auxiliary variable in flow package +description keyword to specify the name of an auxiliary variable in the corresponding flow package. If specified, then the simulated temperatures from this advanced transport package will be copied into the auxiliary variable specified with this name. Note that the flow package must have an auxiliary variable with this name or the program will terminate with an error. If the flows for this advanced transport package are read from a file, then this option will have no effect. + +block options +name boundnames +type keyword +shape +reader urword +optional true +longname +description REPLACE boundnames {'{#1}': 'well'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'well'} + +block options +name print_temperature +type keyword +reader urword +optional true +longname print calculated temperatures to listing file +description REPLACE print_temperature {'{#1}': 'well', '{#2}': 'temperature', '{#3}': 'TEMPERATURE'} + +block options +name print_flows +type keyword +reader urword +optional true +longname print calculated flows to listing file +description REPLACE print_flows {'{#1}': 'well'} + +block options +name save_flows +type keyword +reader urword +optional true +longname save well flows to budget file +description REPLACE save_flows {'{#1}': 'well'} + +block options +name temperature_filerecord +type record temperature fileout tempfile +shape +reader urword +tagged true +optional true +longname +description + +block options +name temperature +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname stage keyword +description keyword to specify that record corresponds to temperature. + +block options +name tempfile +type string +preserve_case true +shape +in_record true +reader urword +tagged false +optional false +longname file keyword +description name of the binary output file to write temperature information. + +block options +name budget_filerecord +type record budget fileout budgetfile +shape +reader urword +tagged true +optional true +longname +description + +block options +name budget +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname budget keyword +description keyword to specify that record corresponds to the budget. + +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 budgetfile +type string +preserve_case true +shape +in_record true +reader urword +tagged false +optional false +longname file keyword +description name of the binary output file to write budget information. + +block options +name budgetcsv_filerecord +type record budgetcsv fileout budgetcsvfile +shape +reader urword +tagged true +optional true +longname +description + +block options +name budgetcsv +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname budget keyword +description keyword to specify that record corresponds to the budget CSV. + +block options +name budgetcsvfile +type string +preserve_case true +shape +in_record true +reader urword +tagged false +optional false +longname file keyword +description name of the comma-separated value (CSV) output file to write budget summary information. A budget summary record will be written to this file for each time step of the simulation. + +block options +name ts_filerecord +type record ts6 filein ts6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name ts6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname head keyword +description keyword to specify that record corresponds to a time-series file. + +block options +name filein +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname file keyword +description keyword to specify that an input filename is expected next. + +block options +name ts6_filename +type string +preserve_case true +in_record true +reader urword +optional false +tagged false +longname file name of time series information +description REPLACE timeseriesfile {} + +block options +name obs_filerecord +type record obs6 filein obs6_filename +shape +reader urword +tagged true +optional true +longname +description + +block options +name obs6 +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname obs keyword +description keyword to specify that record corresponds to an observations file. + +block options +name obs6_filename +type string +preserve_case true +in_record true +tagged false +reader urword +optional false +longname obs6 input filename +description REPLACE obs6_filename {'{#1}': 'MWE'} + + +# --------------------- gwe mwe packagedata --------------------- + +block packagedata +name packagedata +type recarray mawno strt ktf fthk aux boundname +shape (maxbound) +reader urword +longname +description + +block packagedata +name mawno +type integer +shape +tagged false +in_record true +reader urword +longname well number for this entry +description integer value that defines the well number associated with the specified PACKAGEDATA data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. Well information must be specified for every well or the program will terminate with an error. The program will also terminate with an error if information for a well is specified more than once. +numeric_index true + +block packagedata +name strt +type double precision +shape +tagged false +in_record true +reader urword +longname starting well temperature +description real value that defines the starting temperature for the well. + +block packagedata +name ktf +type double precision +shape +tagged false +in_record true +reader urword +longname thermal conductivity of the feature +description is the thermal conductivity of the of the interface between the aquifer cell and the feature. + +block packagedata +name fthk +type double precision +shape +tagged false +in_record true +reader urword +longname thickness of the well feature +description real value that defines the thickness of the material through which conduction occurs. Must be greater than 0. + +block packagedata +name aux +type double precision +in_record true +tagged false +shape (naux) +reader urword +time_series true +optional true +longname auxiliary variables +description REPLACE aux {'{#1}': 'well'} + +block packagedata +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname well name +description REPLACE boundname {'{#1}': 'well'} + + +# --------------------- gwe mwe period --------------------- + +block period +name iper +type integer +block_variable True +in_record true +tagged false +shape +valid +reader urword +optional false +longname stress period number +description REPLACE iper {} + +block period +name mweperioddata +type recarray mawno mwesetting +shape +reader urword +longname +description + +block period +name mawno +type integer +shape +tagged false +in_record true +reader urword +longname well number for this entry +description integer value that defines the well number associated with the specified PERIOD data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. +numeric_index true + +block period +name mwesetting +type keystring status temperature rate auxiliaryrecord +shape +tagged false +in_record true +reader urword +longname +description line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWESETTING string include: STATUS, TEMPERATURE, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the temperature of associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated temperature of the well. + +block period +name status +type string +shape +tagged true +in_record true +reader urword +longname well temperature status +description keyword option to define well status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the well. If a well is inactive, then there will be no solute mass fluxes into or out of the well and the inactive value will be written for the well temperature. If a well is constant, then the temperature for the well will be fixed at the user specified value. + +block period +name temperature +type string +shape +tagged true +in_record true +time_series true +reader urword +longname well temperature +description real or character value that defines the temperature for the well. The specified TEMPERATURE is only applied if the well is a constant temperature well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name rate +type string +shape +tagged true +in_record true +reader urword +time_series true +longname well injection temperature +description real or character value that defines the injection solute temperature $^{\circ}C$ for the well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. + +block period +name auxiliaryrecord +type record auxiliary auxname auxval +shape +tagged +in_record true +reader urword +longname +description + +block period +name auxiliary +type keyword +shape +in_record true +reader urword +longname +description keyword for specifying auxiliary variable. + +block period +name auxname +type string +shape +tagged false +in_record true +reader urword +longname +description name for the auxiliary variable to be assigned AUXVAL. AUXNAME must match one of the auxiliary variable names defined in the OPTIONS block. If AUXNAME does not match one of the auxiliary variable names defined in the OPTIONS block the data are ignored. + +block period +name auxval +type double precision +shape +tagged false +in_record true +reader urword +time_series true +longname auxiliary variable value +description value for the auxiliary variable. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. diff --git a/doc/mf6io/mf6ivar/examples/gwe-mwe-example-obs.dat b/doc/mf6io/mf6ivar/examples/gwe-mwe-example-obs.dat new file mode 100644 index 00000000000..a3a8ced3ad9 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwe-mwe-example-obs.dat @@ -0,0 +1,43 @@ +BEGIN options + DIGITS 12 + PRINT_INPUT +END options + +BEGIN continuous FILEOUT gwe_mwe_02.mwe.obs.csv + mwe1mwe MWE 1 1 + mwe2mwe MWE 2 1 + mwe3mwe MWE 3 1 + mwe4mwe MWE 4 1 + mwe1temp TEMPERATURE 1 + mwe2temp TEMPERATURE 2 + mwe3temp TEMPERATURE 3 + mwe4temp TEMPERATURE 4 + mwe1stor STORAGE 1 + mwe2stor STORAGE 2 + mwe3stor STORAGE 3 + mwe4stor STORAGE 4 + mwe1cnst CONSTANT 1 + mwe2cnst CONSTANT 2 + mwe3cnst CONSTANT 3 + mwe4cnst CONSTANT 4 + mwe1fmvr FROM-MVR 1 + mwe2fmvr FROM-MVR 2 + mwe3fmvr FROM-MVR 3 + mwe4fmvr FROM-MVR 4 + mwe1rate RATE 1 + mwe2rate RATE 2 + mwe3rate RATE 3 + mwe4rate RATE 4 + mwe1rtmv RATE-TO-MVR 1 + mwe2rtmv RATE-TO-MVR 2 + mwe3rtmv RATE-TO-MVR 3 + mwe4rtmv RATE-TO-MVR 4 + mwe1fwrt FW-RATE 1 + mwe2fwrt FW-RATE 2 + mwe3fwrt FW-RATE 3 + mwe4fwrt FW-RATE 4 + mwe1frtm FW-RATE-TO-MVR 1 + mwe2frtm FW-RATE-TO-MVR 2 + mwe3frtm FW-RATE-TO-MVR 3 + mwe4frtm FW-RATE-TO-MVR 4 +END continuous FILEOUT gwe_mwe_02.mwe.obs.csv \ No newline at end of file diff --git a/doc/mf6io/mf6ivar/examples/gwe-mwe-example.dat b/doc/mf6io/mf6ivar/examples/gwe-mwe-example.dat new file mode 100644 index 00000000000..4616f3fa7bc --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwe-mwe-example.dat @@ -0,0 +1,24 @@ +BEGIN OPTIONS + AUXILIARY aux1 aux2 + BOUNDNAMES + PRINT_INPUT + PRINT_TEMPERATURE + PRINT_FLOWS + SAVE_FLOWS + TEMPERATURE FILEOUT gwe_mwe_02.mwe.bin + BUDGET FILEOUT gwe_mwe_02.mwe.bud + OBS6 FILEIN gwe_mwe_02.mwe.obs +END OPTIONS + +BEGIN PACKAGEDATA +# L STRT aux1 aux2 bname + 1 10.00 99.00 999.00 MYWELL1 + 2 10.00 99.00 999.00 MYWELL2 + 3 10.00 99.00 999.00 MYWELL3 +END PACKAGEDATA + +BEGIN PERIOD 1 + 1 STATUS ACTIVE + 2 STATUS ACTIVE + 3 STATUS ACTIVE +END PERIOD 1 diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index e68ea78a611..e30c47cf70b 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -1336,6 +1336,41 @@ | GWE | LKE | PERIOD | AUXILIARY | KEYWORD | keyword for specifying auxiliary variable. | | GWE | LKE | PERIOD | AUXNAME | STRING | name for the auxiliary variable to be assigned AUXVAL. AUXNAME must match one of the auxiliary variable names defined in the OPTIONS block. If AUXNAME does not match one of the auxiliary variable names defined in the OPTIONS block the data are ignored. | | GWE | LKE | PERIOD | AUXVAL | DOUBLE PRECISION | value for the auxiliary variable. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | +| GWE | MWE | OPTIONS | FLOW_PACKAGE_NAME | STRING | keyword to specify the name of the corresponding flow package. If not specified, then the corresponding flow package must have the same name as this advanced transport package (the name associated with this package in the GWE name file). | +| GWE | MWE | OPTIONS | AUXILIARY | STRING (NAUX) | defines an array of one or more auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided on this line; however, lists of information provided in subsequent blocks must have a column of data for each auxiliary variable name defined here. The number of auxiliary variables detected on this line determines the value for naux. Comments cannot be provided anywhere on this line as they will be interpreted as auxiliary variable names. Auxiliary variables may not be used by the package, but they will be available for use by other parts of the program. The program will terminate with an error if auxiliary variables are specified on more than one line in the options block. | +| GWE | MWE | OPTIONS | FLOW_PACKAGE_AUXILIARY_NAME | STRING | keyword to specify the name of an auxiliary variable in the corresponding flow package. If specified, then the simulated temperatures from this advanced transport package will be copied into the auxiliary variable specified with this name. Note that the flow package must have an auxiliary variable with this name or the program will terminate with an error. If the flows for this advanced transport package are read from a file, then this option will have no effect. | +| GWE | MWE | OPTIONS | BOUNDNAMES | KEYWORD | keyword to indicate that boundary names may be provided with the list of well cells. | +| GWE | MWE | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of well information will be written to the listing file immediately after it is read. | +| GWE | MWE | OPTIONS | PRINT_TEMPERATURE | KEYWORD | keyword to indicate that the list of well temperature will be printed to the listing file for every stress period in which ``TEMPERATURE PRINT'' is specified in Output Control. If there is no Output Control option and PRINT\_TEMPERATURE is specified, then temperature are printed for the last time step of each stress period. | +| GWE | MWE | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of well flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. | +| GWE | MWE | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that well flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | +| GWE | MWE | OPTIONS | TEMPERATURE | KEYWORD | keyword to specify that record corresponds to temperature. | +| GWE | MWE | OPTIONS | TEMPFILE | STRING | name of the binary output file to write temperature information. | +| GWE | MWE | OPTIONS | BUDGET | KEYWORD | keyword to specify that record corresponds to the budget. | +| GWE | MWE | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. | +| GWE | MWE | OPTIONS | BUDGETFILE | STRING | name of the binary output file to write budget information. | +| GWE | MWE | OPTIONS | BUDGETCSV | KEYWORD | keyword to specify that record corresponds to the budget CSV. | +| GWE | MWE | OPTIONS | BUDGETCSVFILE | STRING | name of the comma-separated value (CSV) output file to write budget summary information. A budget summary record will be written to this file for each time step of the simulation. | +| GWE | MWE | OPTIONS | TS6 | KEYWORD | keyword to specify that record corresponds to a time-series file. | +| GWE | MWE | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | +| GWE | MWE | OPTIONS | TS6_FILENAME | STRING | defines a time-series file defining time series that can be used to assign time-varying values. See the ``Time-Variable Input'' section for instructions on using the time-series capability. | +| GWE | MWE | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. | +| GWE | MWE | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the MWE package. See the ``Observation utility'' section for instructions for preparing observation input files. Tables \ref{table:gwf-obstypetable} and \ref{table:gwt-obstypetable} lists observation type(s) supported by the MWE package. | +| GWE | MWE | PACKAGEDATA | MAWNO | INTEGER | integer value that defines the well number associated with the specified PACKAGEDATA data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. Well information must be specified for every well or the program will terminate with an error. The program will also terminate with an error if information for a well is specified more than once. | +| GWE | MWE | PACKAGEDATA | STRT | DOUBLE PRECISION | real value that defines the starting temperature for the well. | +| GWE | MWE | PACKAGEDATA | KTF | DOUBLE PRECISION | is the thermal conductivity of the of the interface between the aquifer cell and the feature. | +| GWE | MWE | PACKAGEDATA | FTHK | DOUBLE PRECISION | real value that defines the thickness of the material through which conduction occurs. Must be greater than 0. | +| GWE | MWE | PACKAGEDATA | AUX | DOUBLE PRECISION (NAUX) | represents the values of the auxiliary variables for each well. The values of auxiliary variables must be present for each well. The values must be specified in the order of the auxiliary variables specified in the OPTIONS block. If the package supports time series and the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | +| GWE | MWE | PACKAGEDATA | BOUNDNAME | STRING | name of the well cell. BOUNDNAME is an ASCII character variable that can contain as many as 40 characters. If BOUNDNAME contains spaces in it, then the entire name must be enclosed within single quotes. | +| GWE | MWE | PERIOD | IPER | INTEGER | integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. | +| GWE | MWE | PERIOD | MAWNO | INTEGER | integer value that defines the well number associated with the specified PERIOD data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. | +| GWE | MWE | PERIOD | MWESETTING | KEYSTRING | line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWESETTING string include: STATUS, TEMPERATURE, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the temperature of associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated temperature of the well. | +| GWE | MWE | PERIOD | STATUS | STRING | keyword option to define well status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the well. If a well is inactive, then there will be no solute mass fluxes into or out of the well and the inactive value will be written for the well temperature. If a well is constant, then the temperature for the well will be fixed at the user specified value. | +| GWE | MWE | PERIOD | TEMPERATURE | STRING | real or character value that defines the temperature for the well. The specified TEMPERATURE is only applied if the well is a constant temperature well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | +| GWE | MWE | PERIOD | RATE | STRING | real or character value that defines the injection solute temperature $^{\circ}C$ for the well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | +| GWE | MWE | PERIOD | AUXILIARY | KEYWORD | keyword for specifying auxiliary variable. | +| GWE | MWE | PERIOD | AUXNAME | STRING | name for the auxiliary variable to be assigned AUXVAL. AUXNAME must match one of the auxiliary variable names defined in the OPTIONS block. If AUXNAME does not match one of the auxiliary variable names defined in the OPTIONS block the data are ignored. | +| GWE | MWE | PERIOD | AUXVAL | DOUBLE PRECISION | value for the auxiliary variable. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | | GWE | NAM | OPTIONS | LIST | STRING | is name of the listing file to create for this GWE model. If not specified, then the name of the list file will be the basename of the GWE model name file and the ``.lst'' extension. For example, if the GWE name file is called ``my.model.nam'' then the list file will be called ``my.model.lst''. | | GWE | NAM | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of all model stress package information will be written to the listing file immediately after it is read. | | GWE | NAM | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of all model package flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. | diff --git a/doc/mf6io/mf6ivar/mf6ivar.py b/doc/mf6io/mf6ivar/mf6ivar.py index 25d59480ac5..1170f4a7428 100644 --- a/doc/mf6io/mf6ivar/mf6ivar.py +++ b/doc/mf6io/mf6ivar/mf6ivar.py @@ -749,6 +749,7 @@ def write_appendix(texdir, allblocks): "gwe-est", "gwe-ic", "gwe-lke", + "gwe-mwe", "gwe-nam", "gwe-oc", "gwe-ssm", diff --git a/doc/mf6io/mf6ivar/tex/appendixA.tex b/doc/mf6io/mf6ivar/tex/appendixA.tex index ccdf2aab7a3..cf214f3156f 100644 --- a/doc/mf6io/mf6ivar/tex/appendixA.tex +++ b/doc/mf6io/mf6ivar/tex/appendixA.tex @@ -284,6 +284,10 @@ GWE & LKE & PACKAGEDATA & yes \\ GWE & LKE & PERIOD & yes \\ \hline +GWE & MWE & OPTIONS & yes \\ +GWE & MWE & PACKAGEDATA & yes \\ +GWE & MWE & PERIOD & yes \\ +\hline GWE & NAM & OPTIONS & yes \\ GWE & NAM & PACKAGES & yes \\ \hline diff --git a/doc/mf6io/mf6ivar/tex/gwe-mwe-desc.tex b/doc/mf6io/mf6ivar/tex/gwe-mwe-desc.tex new file mode 100644 index 00000000000..04ee7036eb2 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mwe-desc.tex @@ -0,0 +1,92 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\item \texttt{flow\_package\_name}---keyword to specify the name of the corresponding flow package. If not specified, then the corresponding flow package must have the same name as this advanced transport package (the name associated with this package in the GWE name file). + +\item \texttt{auxiliary}---defines an array of one or more auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided on this line; however, lists of information provided in subsequent blocks must have a column of data for each auxiliary variable name defined here. The number of auxiliary variables detected on this line determines the value for naux. Comments cannot be provided anywhere on this line as they will be interpreted as auxiliary variable names. Auxiliary variables may not be used by the package, but they will be available for use by other parts of the program. The program will terminate with an error if auxiliary variables are specified on more than one line in the options block. + +\item \texttt{flow\_package\_auxiliary\_name}---keyword to specify the name of an auxiliary variable in the corresponding flow package. If specified, then the simulated temperatures from this advanced transport package will be copied into the auxiliary variable specified with this name. Note that the flow package must have an auxiliary variable with this name or the program will terminate with an error. If the flows for this advanced transport package are read from a file, then this option will have no effect. + +\item \texttt{BOUNDNAMES}---keyword to indicate that boundary names may be provided with the list of well cells. + +\item \texttt{PRINT\_INPUT}---keyword to indicate that the list of well information will be written to the listing file immediately after it is read. + +\item \texttt{PRINT\_TEMPERATURE}---keyword to indicate that the list of well temperature will be printed to the listing file for every stress period in which ``TEMPERATURE PRINT'' is specified in Output Control. If there is no Output Control option and PRINT\_TEMPERATURE is specified, then temperature are printed for the last time step of each stress period. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that the list of well flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. + +\item \texttt{SAVE\_FLOWS}---keyword to indicate that well flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. + +\item \texttt{TEMPERATURE}---keyword to specify that record corresponds to temperature. + +\item \texttt{tempfile}---name of the binary output file to write temperature information. + +\item \texttt{BUDGET}---keyword to specify that record corresponds to the budget. + +\item \texttt{FILEOUT}---keyword to specify that an output filename is expected next. + +\item \texttt{budgetfile}---name of the binary output file to write budget information. + +\item \texttt{BUDGETCSV}---keyword to specify that record corresponds to the budget CSV. + +\item \texttt{budgetcsvfile}---name of the comma-separated value (CSV) output file to write budget summary information. A budget summary record will be written to this file for each time step of the simulation. + +\item \texttt{TS6}---keyword to specify that record corresponds to a time-series file. + +\item \texttt{FILEIN}---keyword to specify that an input filename is expected next. + +\item \texttt{ts6\_filename}---defines a time-series file defining time series that can be used to assign time-varying values. See the ``Time-Variable Input'' section for instructions on using the time-series capability. + +\item \texttt{OBS6}---keyword to specify that record corresponds to an observations file. + +\item \texttt{obs6\_filename}---name of input file to define observations for the MWE package. See the ``Observation utility'' section for instructions for preparing observation input files. Tables \ref{table:gwf-obstypetable} and \ref{table:gwt-obstypetable} lists observation type(s) supported by the MWE package. + +\end{description} +\item \textbf{Block: PACKAGEDATA} + +\begin{description} +\item \texttt{mawno}---integer value that defines the well number associated with the specified PACKAGEDATA data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. Well information must be specified for every well or the program will terminate with an error. The program will also terminate with an error if information for a well is specified more than once. + +\item \texttt{strt}---real value that defines the starting temperature for the well. + +\item \texttt{ktf}---is the thermal conductivity of the of the interface between the aquifer cell and the feature. + +\item \texttt{fthk}---real value that defines the thickness of the material through which conduction occurs. Must be greater than 0. + +\item \textcolor{blue}{\texttt{aux}---represents the values of the auxiliary variables for each well. The values of auxiliary variables must be present for each well. The values must be specified in the order of the auxiliary variables specified in the OPTIONS block. If the package supports time series and the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.} + +\item \texttt{boundname}---name of the well cell. BOUNDNAME is an ASCII character variable that can contain as many as 40 characters. If BOUNDNAME contains spaces in it, then the entire name must be enclosed within single quotes. + +\end{description} +\item \textbf{Block: PERIOD} + +\begin{description} +\item \texttt{iper}---integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. + +\item \texttt{mawno}---integer value that defines the well number associated with the specified PERIOD data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. + +\item \texttt{mwesetting}---line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWESETTING string include: STATUS, TEMPERATURE, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the temperature of associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated temperature of the well. + +\begin{lstlisting}[style=blockdefinition] +STATUS +TEMPERATURE <@temperature@> +RATE <@rate@> +AUXILIARY <@auxval@> +\end{lstlisting} + +\item \texttt{status}---keyword option to define well status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the well. If a well is inactive, then there will be no solute mass fluxes into or out of the well and the inactive value will be written for the well temperature. If a well is constant, then the temperature for the well will be fixed at the user specified value. + +\item \textcolor{blue}{\texttt{temperature}---real or character value that defines the temperature for the well. The specified TEMPERATURE is only applied if the well is a constant temperature well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.} + +\item \textcolor{blue}{\texttt{rate}---real or character value that defines the injection solute temperature $^{\circ}C$ for the well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.} + +\item \texttt{AUXILIARY}---keyword for specifying auxiliary variable. + +\item \texttt{auxname}---name for the auxiliary variable to be assigned AUXVAL. AUXNAME must match one of the auxiliary variable names defined in the OPTIONS block. If AUXNAME does not match one of the auxiliary variable names defined in the OPTIONS block the data are ignored. + +\item \textcolor{blue}{\texttt{auxval}---value for the auxiliary variable. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.} + +\end{description} + diff --git a/doc/mf6io/mf6ivar/tex/gwe-mwe-options.dat b/doc/mf6io/mf6ivar/tex/gwe-mwe-options.dat new file mode 100644 index 00000000000..ef83c7fa718 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mwe-options.dat @@ -0,0 +1,15 @@ +BEGIN OPTIONS + [FLOW_PACKAGE_NAME ] + [AUXILIARY ] + [FLOW_PACKAGE_AUXILIARY_NAME ] + [BOUNDNAMES] + [PRINT_INPUT] + [PRINT_TEMPERATURE] + [PRINT_FLOWS] + [SAVE_FLOWS] + [TEMPERATURE FILEOUT ] + [BUDGET FILEOUT ] + [BUDGETCSV FILEOUT ] + [TS6 FILEIN ] + [OBS6 FILEIN ] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/gwe-mwe-packagedata.dat b/doc/mf6io/mf6ivar/tex/gwe-mwe-packagedata.dat new file mode 100644 index 00000000000..0d46cfac181 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mwe-packagedata.dat @@ -0,0 +1,5 @@ +BEGIN PACKAGEDATA + [<@aux(naux)@>] [] + [<@aux(naux)@>] [] + ... +END PACKAGEDATA diff --git a/doc/mf6io/mf6ivar/tex/gwe-mwe-period.dat b/doc/mf6io/mf6ivar/tex/gwe-mwe-period.dat new file mode 100644 index 00000000000..4c9a183703a --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mwe-period.dat @@ -0,0 +1,5 @@ +BEGIN PERIOD + + + ... +END PERIOD diff --git a/doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex b/doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex index a3230eee906..1f558f3244e 100644 --- a/doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex +++ b/doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex @@ -35,7 +35,7 @@ \item \texttt{irch}---IRCH is the layer number that defines the layer in each vertical column where recharge is applied. If IRCH is omitted, recharge by default is applied to cells in layer 1. IRCH can only be used if READASARRAYS is specified in the OPTIONS block. If IRCH is specified, it must be specified as the first variable in the PERIOD block or MODFLOW will terminate with an error. -\item \textcolor{blue}{\texttt{recharge}---is the recharge flux rate ($LT^{-1}$). This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the "Using Time-Array Series in a Package" section).} +\item \textcolor{blue}{\texttt{recharge}---is the recharge flux rate ($LT^{-1}$). This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the ``Using Time-Array Series in a Package'' section).} \item \textcolor{blue}{\texttt{aux}---is an array of values for auxiliary variable aux(iaux), where iaux is a value from 1 to naux, and aux(iaux) must be listed as part of the auxiliary variables. A separate array can be specified for each auxiliary variable. If an array is not specified for an auxiliary variable, then a value of zero is assigned. If the value specified here for the auxiliary variable is the same as auxmultname, then the recharge array will be multiplied by this array.} diff --git a/make/makefile b/make/makefile index 068f81bd825..07f580b3900 100644 --- a/make/makefile +++ b/make/makefile @@ -280,6 +280,7 @@ $(OBJDIR)/GhostNode.o \ $(OBJDIR)/gwf3evt8.o \ $(OBJDIR)/gwf3chd8.o \ $(OBJDIR)/gwe1sfe1.o \ +$(OBJDIR)/gwe1mwe1.o \ $(OBJDIR)/gwe1lke1.o \ $(OBJDIR)/gwe1est1.o \ $(OBJDIR)/gwe1esl1.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 814dc568652..ad8ea63a9d0 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -153,6 +153,7 @@ + diff --git a/src/Model/GroundWaterEnergy/gwe1.f90 b/src/Model/GroundWaterEnergy/gwe1.f90 index e2c86d2c738..dd5536dba78 100644 --- a/src/Model/GroundWaterEnergy/gwe1.f90 +++ b/src/Model/GroundWaterEnergy/gwe1.f90 @@ -764,7 +764,7 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & use GweEslModule, only: esl_create use GweLkeModule, only: lke_create use GweSfeModule, only: sfe_create - !use GweMweModule, only: mwe_create + use GweMweModule, only: mwe_create !use GweUzeModule, only: uze_create use ApiModule, only: api_create ! -- dummy @@ -798,10 +798,10 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & call sfe_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & pakname, this%fmi, this%eqnsclfac, this%gwecommon, & this%depvartype, this%depvarunit, this%depvarunitabbrev) - !case ('MWE6') - ! call mwe_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & - ! pakname, this%fmi, this%tsplab, this%eqnsclfac, & - ! this%gwecommon) + case ('MWE6') + call mwe_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & + pakname, this%fmi, this%eqnsclfac, this%gwecommon, & + this%depvartype, this%depvarunit, this%depvarunitabbrev) !case ('UZE6') ! call uze_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & ! pakname, this%fmi, this%tsplab, this%eqnsclfac, & diff --git a/src/Model/GroundWaterEnergy/gwe1mwe1.f90 b/src/Model/GroundWaterEnergy/gwe1mwe1.f90 new file mode 100644 index 00000000000..4b6bb76a38d --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe1mwe1.f90 @@ -0,0 +1,1018 @@ +! -- Multi-Aquifer Well Energy Transport Module +! -- todo: save the mwe temperature into the mwe aux variable? +! -- todo: calculate the maw DENSE aux variable using temperature? +! +! MAW flows (flowbudptr) index var MWE term Transport Type +!--------------------------------------------------------------------------------- + +! -- terms from MAW that will be handled by parent APT Package +! FLOW-JA-FACE idxbudfjf FLOW-JA-FACE cv2cv (note that this doesn't exist for MAW) +! GWF (aux FLOW-AREA) idxbudgwf GWF cv2gwf +! STORAGE (aux VOLUME) idxbudsto none used for cv volumes +! FROM-MVR idxbudfmvr FROM-MVR q * text = this%qfrommvr(:) +! TO-MVR idxbudtmvr TO-MVR q * tfeat + +! -- MAW terms +! RATE idxbudrate RATE q < 0: q * twell, else q * tuser +! FW-RATE idxbudfwrt FW-RATE q * twell +! RATE-TO-MVR idxbudrtmv RATE-TO-MVR q * twell +! FW-RATE-TO-MVR idxbudfrtm FW-RATE-TO-MVR q * twell +! WELL-AQUIFER CONDUCTION idxbudmwcd MW-CONDUCTION K_t_f * WetArea / thickness + +! -- terms from MAW that should be skipped +! CONSTANT-TO-MVR ? CONSTANT-TO-MVR q * twell + +! -- terms from a flow file that should be skipped +! CONSTANT none none none +! AUXILIARY none none none + +! -- terms that are written to the energy transport budget file +! none none STORAGE (aux ENER) dE/dt +! none none AUXILIARY none +! none none CONSTANT accumulate +! +! +module GweMweModule + + use KindModule, only: DP, I4B + use ConstantsModule, only: DZERO, LINELENGTH + use SimModule, only: store_error + use BndModule, only: BndType, GetBndFromList + use TspFmiModule, only: TspFmiType + use MawModule, only: MawType + use ObserveModule, only: ObserveType + use TspAptModule, only: TspAptType, apt_process_obsID, & + apt_process_obsID12 + use GweInputDataModule, only: GweInputDataType + use MatrixBaseModule + + implicit none + + public mwe_create + + character(len=*), parameter :: ftype = 'MWE' + character(len=*), parameter :: flowtype = 'MAW' + character(len=16) :: text = ' MWE' + + type, extends(TspAptType) :: GweMweType + + type(GweInputDataType), pointer :: gwecommon => null() !< pointer to shared gwe data used by multiple packages but set in mst + + integer(I4B), pointer :: idxbudrate => null() ! index of well rate terms in flowbudptr + integer(I4B), pointer :: idxbudfwrt => null() ! index of flowing well rate terms in flowbudptr + integer(I4B), pointer :: idxbudrtmv => null() ! index of rate to mover terms in flowbudptr + integer(I4B), pointer :: idxbudfrtm => null() ! index of flowing well rate to mover terms in flowbudptr + integer(I4B), pointer :: idxbudmwcd => null() ! index of well bore conduction terms in flowbudptr + real(DP), dimension(:), pointer, contiguous :: temprate => null() ! well rate temperature + + contains + + procedure :: bnd_da => mwe_da + procedure :: allocate_scalars + procedure :: apt_allocate_arrays => mwe_allocate_arrays + procedure :: find_apt_package => find_mwe_package + procedure :: pak_fc_expanded => mwe_fc_expanded + procedure :: pak_solve => mwe_solve + procedure :: pak_get_nbudterms => mwe_get_nbudterms + procedure :: pak_setup_budobj => mwe_setup_budobj + procedure :: pak_fill_budobj => mwe_fill_budobj + procedure :: mwe_rate_term + procedure :: mwe_fwrt_term + procedure :: mwe_rtmv_term + procedure :: mwe_frtm_term + procedure :: pak_df_obs => mwe_df_obs + procedure :: pak_rp_obs => mwe_rp_obs + procedure :: pak_bd_obs => mwe_bd_obs + procedure :: pak_set_stressperiod => mwe_set_stressperiod + + end type GweMweType + +contains + + !> Create new MWE package + !< + subroutine mwe_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname, & + fmi, eqnsclfac, gwecommon, dvt, dvu, dvua) + ! -- dummy + class(BndType), pointer :: packobj + integer(I4B), intent(in) :: id + integer(I4B), intent(in) :: ibcnum + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + character(len=*), intent(in) :: namemodel + character(len=*), intent(in) :: pakname + type(TspFmiType), pointer :: fmi + real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor + type(GweInputDataType), intent(in), target :: gwecommon !< shared data container for use by multiple GWE packages + character(len=*), intent(in) :: dvt !< For GWE, set to "TEMPERATURE" in TspAptType + character(len=*), intent(in) :: dvu !< For GWE, set to "energy" in TspAptType + character(len=*), intent(in) :: dvua !< For GWE, set to "E" in TspAptType + ! -- local + type(GweMweType), pointer :: mweobj + ! + ! -- Allocate the object and assign values to object variables + allocate (mweobj) + packobj => mweobj + ! + ! -- Create name and memory path + call packobj%set_names(ibcnum, namemodel, pakname, ftype) + packobj%text = text + ! + ! -- Allocate scalars + call mweobj%allocate_scalars() + ! + ! -- Initialize package + call packobj%pack_initialize() + ! + packobj%inunit = inunit + packobj%iout = iout + packobj%id = id + packobj%ibcnum = ibcnum + packobj%ncolbnd = 1 + packobj%iscloc = 1 + ! + ! -- Store pointer to flow model interface. When the GwfGwe exchange is + ! created, it sets fmi%bndlist so that the GWE model has access to all + ! the flow packages + mweobj%fmi => fmi + ! + ! -- Store pointer to governing equation scale factor + mweobj%eqnsclfac => eqnsclfac + ! + ! -- Store pointer to shared data module for accessing cpw, rhow + ! for the budget calculations, and for accessing the latent heat of + ! vaporization for evaporative cooling. + mweobj%gwecommon => gwecommon + ! + ! -- Set labels that will be used in generalized APT class + mweobj%depvartype = dvt + mweobj%depvarunit = dvu + mweobj%depvarunitabbrev = dvua + ! + ! -- Return + return + end subroutine mwe_create + + !> @brief Find corresponding mwe package + !< + subroutine find_mwe_package(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweMweType) :: this + ! -- local + character(len=LINELENGTH) :: errmsg + class(BndType), pointer :: packobj + integer(I4B) :: ip, icount + integer(I4B) :: nbudterm + logical :: found + ! + ! -- Initialize found to false, and error later if flow package cannot + ! be found + found = .false. + ! + ! -- If user is specifying flows in a binary budget file, then set up + ! the budget file reader, otherwise set a pointer to the flow package + ! budobj + if (this%fmi%flows_from_file) then + call this%fmi%set_aptbudobj_pointer(this%flowpackagename, this%flowbudptr) + if (associated(this%flowbudptr)) found = .true. + ! + else + if (associated(this%fmi%gwfbndlist)) then + ! -- Look through gwfbndlist for a flow package with the same name as + ! this transport package name + do ip = 1, this%fmi%gwfbndlist%Count() + packobj => GetBndFromList(this%fmi%gwfbndlist, ip) + if (packobj%packName == this%flowpackagename) then + found = .true. + ! + ! -- Store BndType pointer to packobj, and then + ! use the select type to point to the budobj in flow package + this%flowpackagebnd => packobj + select type (packobj) + type is (MawType) + this%flowbudptr => packobj%budobj + end select + end if + if (found) exit + end do + end if + end if + ! + ! -- Error if flow package not found + if (.not. found) then + write (errmsg, '(a)') 'COULD NOT FIND FLOW PACKAGE WITH NAME '& + &//trim(adjustl(this%flowpackagename))//'.' + call store_error(errmsg) + call this%parser%StoreErrorUnit() + end if + ! + ! -- Allocate space for idxbudssm, which indicates whether this is a + ! special budget term or one that is a general source and sink + nbudterm = this%flowbudptr%nbudterm + call mem_allocate(this%idxbudssm, nbudterm, 'IDXBUDSSM', this%memoryPath) + ! + ! -- Process budget terms and identify special budget terms + write (this%iout, '(/, a, a)') & + 'PROCESSING '//ftype//' INFORMATION FOR ', this%packName + write (this%iout, '(a)') ' IDENTIFYING FLOW TERMS IN '//flowtype//' PACKAGE' + write (this%iout, '(a, i0)') & + ' NUMBER OF '//flowtype//' = ', this%flowbudptr%ncv + icount = 1 + do ip = 1, this%flowbudptr%nbudterm + select case (trim(adjustl(this%flowbudptr%budterm(ip)%flowtype))) + case ('FLOW-JA-FACE') + this%idxbudfjf = ip + this%idxbudssm(ip) = 0 + case ('GWF') + this%idxbudgwf = ip + this%idxbudssm(ip) = 0 + case ('STORAGE') + this%idxbudsto = ip + this%idxbudssm(ip) = 0 + case ('RATE') + this%idxbudrate = ip + this%idxbudssm(ip) = 0 + case ('FW-RATE') + this%idxbudfwrt = ip + this%idxbudssm(ip) = 0 + case ('RATE-TO-MVR') + this%idxbudrtmv = ip + this%idxbudssm(ip) = 0 + case ('FW-RATE-TO-MVR') + this%idxbudfrtm = ip + this%idxbudssm(ip) = 0 + case ('TO-MVR') + this%idxbudtmvr = ip + this%idxbudssm(ip) = 0 + case ('FROM-MVR') + this%idxbudfmvr = ip + this%idxbudssm(ip) = 0 + case ('AUXILIARY') + this%idxbudaux = ip + this%idxbudssm(ip) = 0 + case default + ! + ! -- Set idxbudssm equal to a column index for where the tempeartures + ! are stored in the tempbud(nbudssm, ncv) array + this%idxbudssm(ip) = icount + icount = icount + 1 + end select + write (this%iout, '(a, i0, " = ", a,/, a, i0)') & + ' TERM ', ip, trim(adjustl(this%flowbudptr%budterm(ip)%flowtype)), & + ' MAX NO. OF ENTRIES = ', this%flowbudptr%budterm(ip)%maxlist + end do + write (this%iout, '(a, //)') 'DONE PROCESSING '//ftype//' INFORMATION' + ! + ! -- Streambed conduction term + this%idxbudmwcd = this%idxbudgwf + ! + ! -- Return + return + end subroutine find_mwe_package + + !> @brief Add matrix terms related to MWE + !! + !! This routine is called from TspAptType%apt_fc_expanded() in + !! order to add matrix terms specifically for MWE + !< + subroutine mwe_fc_expanded(this, rhs, ia, idxglo, matrix_sln) + ! -- dummy + class(GweMweType) :: this + real(DP), dimension(:), intent(inout) :: rhs + integer(I4B), dimension(:), intent(in) :: ia + integer(I4B), dimension(:), intent(in) :: idxglo + class(MatrixBaseType), pointer :: matrix_sln + ! -- local + integer(I4B) :: j, n, n1, n2 + integer(I4B) :: iloc + integer(I4B) :: iposd, iposoffd + integer(I4B) :: ipossymd, ipossymoffd + integer(I4B) :: auxpos + real(DP) :: rrate + real(DP) :: rhsval + real(DP) :: hcofval + real(DP) :: ctherm ! kluge? + real(DP) :: wa !< wetted area + real(DP) :: ktf !< thermal conductivity of streambed material + real(DP) :: s !< thickness of conductive wellbore material + ! + ! -- Add puping rate contribution + if (this%idxbudrate /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrate)%nlist + call this%mwe_rate_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) + iposd = this%idxpakdiag(n1) + call matrix_sln%add_value_pos(iposd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add flowing well rate contribution + if (this%idxbudfwrt /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudfwrt)%nlist + call this%mwe_fwrt_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) + iposd = this%idxpakdiag(n1) + call matrix_sln%add_value_pos(iposd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add rate to mover contribution + if (this%idxbudrtmv /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrtmv)%nlist + call this%mwe_rtmv_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) + iposd = this%idxpakdiag(n1) + call matrix_sln%add_value_pos(iposd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add puping rate contribution + if (this%idxbudfrtm /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudfrtm)%nlist + call this%mwe_frtm_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) + iposd = this%idxpakdiag(n1) + call matrix_sln%add_value_pos(iposd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add wellbore conduction contribution + do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + ! + ! -- Set n to feature number and process if active features + n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j) + if (this%iboundpak(n) /= 0) then + ! + ! -- Set acoef and rhs to negative so they are relative to mwe and not gwe + auxpos = this%flowbudptr%budterm(this%idxbudgwf)%naux + wa = this%flowbudptr%budterm(this%idxbudgwf)%auxvar(auxpos, j) + ktf = this%ktf(n) + s = this%rfeatthk(n) + ctherm = ktf * wa / s + ! + ! -- Add to mwe row + iposd = this%idxdglo(j) + iposoffd = this%idxoffdglo(j) + call matrix_sln%add_value_pos(iposd, -ctherm) ! kluge note: make sure the signs on ctherm are correct here and below + call matrix_sln%add_value_pos(iposoffd, ctherm) + ! + ! -- Add to gwe row for mwe connection + ipossymd = this%idxsymdglo(j) + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymd, -ctherm) + call matrix_sln%add_value_pos(ipossymoffd, ctherm) + end if + end do + ! + ! -- Return + return + end subroutine mwe_fc_expanded + + !> @brief Add terms specific to multi-aquifer wells to the explicit multi- + !! aquifer well energy transport solve + !< + subroutine mwe_solve(this) + ! -- dummy + class(GweMweType) :: this + ! -- local + integer(I4B) :: j + integer(I4B) :: n1, n2 + real(DP) :: rrate + ! + ! -- Add well pumping contribution + if (this%idxbudrate /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrate)%nlist + call this%mwe_rate_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add flowing well rate contribution + if (this%idxbudfwrt /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudfwrt)%nlist + call this%mwe_fwrt_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add well pumping rate to mover contribution + if (this%idxbudrtmv /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrtmv)%nlist + call this%mwe_rtmv_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add flowing well rate to mover contribution + if (this%idxbudfrtm /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudfrtm)%nlist + call this%mwe_frtm_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Return + return + end subroutine mwe_solve + + !> @brief Function to return the number of budget terms just for this package + !! + !! This overrides a function in the parent class. + !< + function mwe_get_nbudterms(this) result(nbudterms) + ! -- dummy + class(GweMweType) :: this + ! -- return + integer(I4B) :: nbudterms + ! + ! -- Number of potential budget terms is 5 + nbudterms = 1 ! RATE + if (this%idxbudfwrt /= 0) nbudterms = nbudterms + 1 + if (this%idxbudrtmv /= 0) nbudterms = nbudterms + 1 + if (this%idxbudfrtm /= 0) nbudterms = nbudterms + 1 + if (this%idxbudmwcd /= 0) nbudterms = nbudterms + 1 + ! + ! -- Return + return + end function mwe_get_nbudterms + + !> @brief Set up the budget object that stores all the mwe flows + !< + subroutine mwe_setup_budobj(this, idx) + ! -- modules + use ConstantsModule, only: LENBUDTXT + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(inout) :: idx + ! -- local + integer(I4B) :: n, n1, n2 + integer(I4B) :: maxlist, naux + real(DP) :: q + character(len=LENBUDTXT) :: text + ! + ! -- User-specified rate + text = ' RATE' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudrate)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + ! + ! -- Flowing well rate + if (this%idxbudfwrt /= 0) then + text = ' FW-RATE' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudfwrt)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + end if + ! + ! -- User-specified flow rate to mover + if (this%idxbudrtmv /= 0) then + text = ' RATE-TO-MVR' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudrtmv)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + end if + ! + ! -- Fflowing well rate to mover + if (this%idxbudfrtm /= 0) then + text = ' FW-RATE-TO-MVR' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudfrtm)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + end if + ! + ! -- Conduction through wellbore (and/or filter pack) + text = ' WELLBORE-COND' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudmwcd)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + call this%budobj%budterm(idx)%reset(maxlist) + q = DZERO + do n = 1, maxlist + n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(n) + n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(n) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + end do + ! + ! -- Return + return + end subroutine mwe_setup_budobj + + !> @brief Copy flow terms into this%budobj + !< + subroutine mwe_fill_budobj(this, idx, x, flowja, ccratin, ccratout) + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(inout) :: idx + real(DP), dimension(:), intent(in) :: x + real(DP), dimension(:), contiguous, intent(inout) :: flowja + real(DP), intent(inout) :: ccratin + real(DP), intent(inout) :: ccratout + ! -- local + integer(I4B) :: j, n1, n2 + integer(I4B) :: nlist + integer(I4B) :: igwfnode + integer(I4B) :: idiag + integer(I4B) :: auxpos + real(DP) :: q + real(DP) :: ctherm + real(DP) :: wa !< wetted area + real(DP) :: ktf !< thermal conductivity of streambed material + real(DP) :: s !< thickness of conductive streambed materia + ! + ! -- Rate + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudrate)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%mwe_rate_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + end do + ! + ! -- FW-Rate + if (this%idxbudfwrt /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudfwrt)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%mwe_fwrt_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + end do + end if + ! + ! -- Rate-To-MVR + if (this%idxbudrtmv /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudrtmv)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%mwe_rtmv_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + end do + end if + ! + ! -- FW-Rate-To-MVR + if (this%idxbudfrtm /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudfrtm)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%mwe_frtm_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + end do + end if + ! + ! -- Wellbore-Cond + idx = idx + 1 + call this%budobj%budterm(idx)%reset(this%maxbound) + do j = 1, this%flowbudptr%budterm(this%idxbudmwcd)%nlist + q = DZERO + n1 = this%flowbudptr%budterm(this%idxbudmwcd)%id1(j) + if (this%iboundpak(n1) /= 0) then + igwfnode = this%flowbudptr%budterm(this%idxbudmwcd)%id2(j) + auxpos = this%flowbudptr%budterm(this%idxbudgwf)%naux + wa = this%flowbudptr%budterm(this%idxbudgwf)%auxvar(auxpos, j) + ktf = this%ktf(n1) + s = this%rfeatthk(n1) + ctherm = ktf * wa / s + q = ctherm * (x(igwfnode) - this%xnewpak(n1)) + end if + call this%budobj%budterm(idx)%update_term(n1, igwfnode, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + if (this%iboundpak(n1) /= 0) then + ! -- Contribution to gwe cell budget + this%simvals(j) = this%simvals(j) - q + idiag = this%dis%con%ia(igwfnode) + flowja(idiag) = flowja(idiag) - q + end if + end do + ! + ! -- Return + return + end subroutine mwe_fill_budobj + + !> @brief Allocate scalars specific to the multi-aquifer well energy + !! transport (MWE) package. + !< + subroutine allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweMweType) :: this + ! + ! -- Allocate scalars in TspAptType + call this%TspAptType%allocate_scalars() + ! + ! -- Allocate + call mem_allocate(this%idxbudrate, 'IDXBUDRATE', this%memoryPath) + call mem_allocate(this%idxbudfwrt, 'IDXBUDFWRT', this%memoryPath) + call mem_allocate(this%idxbudrtmv, 'IDXBUDRTMV', this%memoryPath) + call mem_allocate(this%idxbudfrtm, 'IDXBUDFRTM', this%memoryPath) + call mem_allocate(this%idxbudmwcd, 'IDXBUDMWCD', this%memoryPath) + ! + ! -- Initialize + this%idxbudrate = 0 + this%idxbudfwrt = 0 + this%idxbudrtmv = 0 + this%idxbudfrtm = 0 + this%idxbudmwcd = 0 + ! + ! -- Return + return + end subroutine allocate_scalars + + !> @brief Allocate arrays specific to the streamflow mass transport (SFT) + !! package + !< + subroutine mwe_allocate_arrays(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweMweType), intent(inout) :: this + ! -- local + integer(I4B) :: n + ! + ! -- Time series + call mem_allocate(this%temprate, this%ncv, 'TEMPRATE', this%memoryPath) + ! + ! -- Call standard TspAptType allocate arrays + call this%TspAptType%apt_allocate_arrays() + ! + ! -- Initialize + do n = 1, this%ncv + this%temprate(n) = DZERO + end do + ! + ! -- Return + return + end subroutine mwe_allocate_arrays + + !> @brief Deallocate memory associated with MWE package + !< + subroutine mwe_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(GweMweType) :: this + ! + ! -- Deallocate scalars + call mem_deallocate(this%idxbudrate) + call mem_deallocate(this%idxbudfwrt) + call mem_deallocate(this%idxbudrtmv) + call mem_deallocate(this%idxbudfrtm) + call mem_deallocate(this%idxbudmwcd) + ! + ! -- Deallocate time series + call mem_deallocate(this%temprate) + ! + ! -- Deallocate scalars in TspAptType + call this%TspAptType%bnd_da() + ! + ! -- Return + return + end subroutine mwe_da + + !> @brief Thermal transport matrix term(s) associcated with a user-specified + !! flow rate (mwe_rate_term) + !< + subroutine mwe_rate_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout), optional :: rrate + real(DP), intent(inout), optional :: rhsval + real(DP), intent(inout), optional :: hcofval + ! -- local + real(DP) :: qbnd + real(DP) :: ctmp + real(DP) :: h, r + ! + n1 = this%flowbudptr%budterm(this%idxbudrate)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudrate)%id2(ientry) + ! -- Note that qbnd is negative for an extracting well + qbnd = this%flowbudptr%budterm(this%idxbudrate)%flow(ientry) + if (qbnd < DZERO) then + ctmp = this%xnewpak(n1) + h = qbnd + r = DZERO + else + ctmp = this%temprate(n1) + h = DZERO + r = -qbnd * ctmp + end if + if (present(rrate)) rrate = qbnd * ctmp * this%eqnsclfac + if (present(rhsval)) rhsval = r * this%eqnsclfac + if (present(hcofval)) hcofval = h * this%eqnsclfac + ! + ! -- Return + return + end subroutine mwe_rate_term + + !> @brief Thermal transport matrix term(s) associcated with a flowing- + !! well rate term associated with pumping (or injection) + !< + subroutine mwe_fwrt_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout), optional :: rrate + real(DP), intent(inout), optional :: rhsval + real(DP), intent(inout), optional :: hcofval + ! -- local + real(DP) :: qbnd + real(DP) :: ctmp + ! + n1 = this%flowbudptr%budterm(this%idxbudfwrt)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudfwrt)%id2(ientry) + qbnd = this%flowbudptr%budterm(this%idxbudfwrt)%flow(ientry) + ctmp = this%xnewpak(n1) + if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac + if (present(rhsval)) rhsval = DZERO + if (present(hcofval)) hcofval = qbnd * this%eqnsclfac + ! + ! -- Return + return + end subroutine mwe_fwrt_term + + !> @brief Thermal transport matrix term(s) associcated with pumped-water- + !! to-mover term (mwe_rtmv_term) + !< + subroutine mwe_rtmv_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout), optional :: rrate + real(DP), intent(inout), optional :: rhsval + real(DP), intent(inout), optional :: hcofval + ! -- local + real(DP) :: qbnd + real(DP) :: ctmp + ! + n1 = this%flowbudptr%budterm(this%idxbudrtmv)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudrtmv)%id2(ientry) + qbnd = this%flowbudptr%budterm(this%idxbudrtmv)%flow(ientry) + ctmp = this%xnewpak(n1) + if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac + if (present(rhsval)) rhsval = DZERO + if (present(hcofval)) hcofval = qbnd * this%eqnsclfac + ! + ! -- Return + return + end subroutine mwe_rtmv_term + + !> @brief Thermal transport matrix term(s) associcated with the flowing- + !! well-rate-to-mover term (mwe_frtm_term) + !< + subroutine mwe_frtm_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweMweType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout), optional :: rrate + real(DP), intent(inout), optional :: rhsval + real(DP), intent(inout), optional :: hcofval + ! -- local + real(DP) :: qbnd + real(DP) :: ctmp + ! + n1 = this%flowbudptr%budterm(this%idxbudfrtm)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudfrtm)%id2(ientry) + qbnd = this%flowbudptr%budterm(this%idxbudfrtm)%flow(ientry) + ctmp = this%xnewpak(n1) + if (present(rrate)) rrate = ctmp * qbnd * this%eqnsclfac + if (present(rhsval)) rhsval = DZERO + if (present(hcofval)) hcofval = qbnd * this%eqnsclfac + ! + ! -- Return + return + end subroutine mwe_frtm_term + + !> @brief Observations + !! + !! Store the observation type supported by the APT package and overide + !! BndType%bnd_df_obs + !< + subroutine mwe_df_obs(this) + ! -- dummy + class(GweMweType) :: this + ! -- local + integer(I4B) :: indx + ! + ! -- Store obs type and assign procedure pointer + ! for temperature observation type. + call this%obs%StoreObsType('temperature', .false., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Flow-ja-face not supported for MWE + !call this%obs%StoreObsType('flow-ja-face', .true., indx) + !this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for from-mvr observation type. + call this%obs%StoreObsType('from-mvr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- To-mvr not supported for mwe + !call this%obs%StoreObsType('to-mvr', .true., indx) + !this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for storage observation type. + call this%obs%StoreObsType('storage', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for constant observation type. + call this%obs%StoreObsType('constant', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type: mwe + call this%obs%StoreObsType('mwe', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID12 + ! + ! -- Store obs type and assign procedure pointer + ! for rate observation type. + call this%obs%StoreObsType('rate', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('fw-rate', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('rate-to-mvr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('fw-rate-to-mvr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Return + return + end subroutine mwe_df_obs + + !> @brief Process package specific obs + !! + !! Method to process specific observations for this package. + !< + subroutine mwe_rp_obs(this, obsrv, found) + ! -- dummy + class(GweMweType), intent(inout) :: this !< package class + type(ObserveType), intent(inout) :: obsrv !< observation object + logical, intent(inout) :: found !< indicate whether observation was found + ! + found = .true. + select case (obsrv%ObsTypeId) + case ('RATE') + call this%rp_obs_byfeature(obsrv) + case ('FW-RATE') + call this%rp_obs_byfeature(obsrv) + case ('RATE-TO-MVR') + call this%rp_obs_byfeature(obsrv) + case ('FW-RATE-TO-MVR') + call this%rp_obs_byfeature(obsrv) + case default + found = .false. + end select + ! + ! -- Return + return + end subroutine mwe_rp_obs + + !> @brief Calculate observation value and pass it back to APT + !< + subroutine mwe_bd_obs(this, obstypeid, jj, v, found) + ! -- dummy + class(GweMweType), intent(inout) :: this + character(len=*), intent(in) :: obstypeid + real(DP), intent(inout) :: v + integer(I4B), intent(in) :: jj + logical, intent(inout) :: found + ! -- local + integer(I4B) :: n1, n2 + ! + found = .true. + select case (obstypeid) + case ('RATE') + if (this%iboundpak(jj) /= 0) then + call this%mwe_rate_term(jj, n1, n2, v) + end if + case ('FW-RATE') + if (this%iboundpak(jj) /= 0 .and. this%idxbudfwrt > 0) then + call this%mwe_fwrt_term(jj, n1, n2, v) + end if + case ('RATE-TO-MVR') + if (this%iboundpak(jj) /= 0 .and. this%idxbudrtmv > 0) then + call this%mwe_rtmv_term(jj, n1, n2, v) + end if + case ('FW-RATE-TO-MVR') + if (this%iboundpak(jj) /= 0 .and. this%idxbudfrtm > 0) then + call this%mwe_frtm_term(jj, n1, n2, v) + end if + case default + found = .false. + end select + ! + ! -- Return + return + end subroutine mwe_bd_obs + + !> @brief Sets the stress period attributes for keyword use. + !< + subroutine mwe_set_stressperiod(this, itemno, keyword, found) + ! -- modules + use TimeSeriesManagerModule, only: read_value_or_time_series_adv + ! -- dummy + class(GweMweType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + character(len=*), intent(in) :: keyword + logical, intent(inout) :: found + ! -- local + character(len=LINELENGTH) :: text + integer(I4B) :: ierr + integer(I4B) :: jj + real(DP), pointer :: bndElem => null() + ! + ! RATE + ! + found = .true. + select case (keyword) + case ('RATE') + ierr = this%apt_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(text) + jj = 1 + bndElem => this%temprate(itemno) + call read_value_or_time_series_adv(text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'RATE') + case default + ! + ! -- Keyword not recognized so return to caller with found = .false. + found = .false. + end select + ! +999 continue + ! + ! -- Return + return + end subroutine mwe_set_stressperiod + +end module GweMweModule diff --git a/src/Model/TransportModel/tsp1apt1.f90 b/src/Model/TransportModel/tsp1apt1.f90 index facc861b3bb..76fc718e033 100644 --- a/src/Model/TransportModel/tsp1apt1.f90 +++ b/src/Model/TransportModel/tsp1apt1.f90 @@ -2813,7 +2813,7 @@ subroutine apt_rp_obs(this) ' must be assigned to a feature with a unique boundname.' call store_error(errmsg) end if - case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE') + case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE') call this%rp_obs_budterm(obsrv, & this%flowbudptr%budterm(this%idxbudgwf)) case ('FLOW-JA-FACE') @@ -2898,7 +2898,7 @@ subroutine apt_bd_obs(this) if (this%iboundpak(jj) /= 0) then v = this%xnewpak(jj) end if - case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE') + case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE') n = this%flowbudptr%budterm(this%idxbudgwf)%id1(jj) if (this%iboundpak(n) /= 0) then igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(jj) diff --git a/src/meson.build b/src/meson.build index 0fd84e272a2..3eca98fdb7e 100644 --- a/src/meson.build +++ b/src/meson.build @@ -71,6 +71,7 @@ modflow_sources = files( 'Model' / 'GroundWaterEnergy' / 'gwe1ic1idm.f90', 'Model' / 'GroundWaterEnergy' / 'gwe1idm.f90', 'Model' / 'GroundWaterEnergy' / 'gwe1lke1.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe1mwe1.f90', 'Model' / 'GroundWaterEnergy' / 'gwe1sfe1.f90', 'Model' / 'GroundWaterFlow' / 'gwf3.f90', 'Model' / 'GroundWaterFlow' / 'gwf3api8.f90',