From 754b59fe2cae3317ed3303f7e8c1297670fb88e9 Mon Sep 17 00:00:00 2001 From: Eric Morway Date: Wed, 17 Apr 2024 10:22:00 -0700 Subject: [PATCH] feat(MVE): adding support for energy mover package --- autotest/test_gwe_mve.py | 872 ++++++++++++++++++ autotest/test_gwe_uze00.py | 8 +- doc/mf6io/gwe/gwe.tex | 4 + doc/mf6io/gwe/mve.tex | 18 + doc/mf6io/mf6ivar/dfn/gwe-mve.dfn | 106 +++ .../mf6ivar/examples/gwe-mve-example.dat | 6 + doc/mf6io/mf6ivar/md/mf6ivar.md | 8 + doc/mf6io/mf6ivar/mf6ivar.py | 1 - doc/mf6io/mf6ivar/tex/appendixA.tex | 2 + doc/mf6io/mf6ivar/tex/gwe-mve-desc.tex | 23 + doc/mf6io/mf6ivar/tex/gwe-mve-options.dat | 7 + src/Exchange/exg-gwegwe.f90 | 2 +- src/Exchange/exg-gwtgwt.f90 | 2 +- src/Model/GroundWaterEnergy/gwe.f90 | 2 +- src/Model/TransportModel/tsp-mvt.f90 | 167 ++-- src/Model/TransportModel/tsp.f90 | 4 +- src/Utilities/BudgetObject.f90 | 12 +- 17 files changed, 1134 insertions(+), 110 deletions(-) create mode 100644 autotest/test_gwe_mve.py create mode 100644 doc/mf6io/gwe/mve.tex create mode 100644 doc/mf6io/mf6ivar/dfn/gwe-mve.dfn create mode 100644 doc/mf6io/mf6ivar/examples/gwe-mve-example.dat create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mve-desc.tex create mode 100644 doc/mf6io/mf6ivar/tex/gwe-mve-options.dat diff --git a/autotest/test_gwe_mve.py b/autotest/test_gwe_mve.py new file mode 100644 index 00000000000..7a6ccdfbdc5 --- /dev/null +++ b/autotest/test_gwe_mve.py @@ -0,0 +1,872 @@ +""" +Uses a sloped surface to test MVR/MVE connections. Connection types include + 1) DRN to SFR + 2) UZF to SFR + 3) DRN to UZF + 4) UZF to UZF + +Stress period 1 shouldn't have any DRN to MVR flows +Stress period 1 should include UZF to UZF and UZF to SFR flows +Stress period 2 shouldn't have any type of TO-MVR flows +Stress period 3 should include DRN to UZF and DRN to SFR flows +Stress period 4 should include all types (listed above) of TO-MVR flows + + (Not to scale: shows model grid design in profile view) + + +-----+ + | +-----+ precip + | | +-----+ runoff | MVR connections: UZF -> SFR + | | | +-----+ ~~~> v UZF -> UZF + | | | | +-----+ ~~~> DRN -> UZF + | | | | | +-----+ ~~~> DRN -> SFR + | | | | | | +-----+ ~~~> + | | | | | | | +-----+ ~~~> SFR channel + | | | | | | | | +-----+ / + | | | | | | | | | gw +-| |-+ + | | | | | | | | |dis Q| +-+ | + +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ + | | | | | | | | | | | + | | | | | | | | | | | + +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ + | | | | | | | | | | | + | | | | | | | | | | | + +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+ + + A note about "gw dis Q" in the 9th column: DRN connections using + MVR/MVE have been set up for all active UZF cells; however, only + the right-most active DRN connection actually transfers groundwater + discharge to land surface to SFR. That is, none of the DRN -> UZF + connections accounting for cascading gw discharge to land surface + should ever be anything but 0.0 (which is verified by the checks) + +""" + +import os + +import flopy +import flopy.utils.binaryfile as bf +import numpy as np +import pytest + +from framework import TestFramework + +include_NWT = False + +cases = ["mve-01"] + +iuz_cell_dict = {} +cell_iuz_dict = {} + +nlay, nrow, ncol = 3, 3, 10 +nper = 4 +perlen = [1.0] * nper +nstp = [1] * nper +tsmult = [1.0] * nper + +delr = 1.0 +delc = 1.0 +strt = 25 + +nouter, ninner = 300, 300 +hclose, rclose, relax = 1e-6, 1e-3, 0.97 + +top = np.stack( + [ + [30.9, 30.8, 30.6, 30.6, 30.5, 30.4, 30.3, 30.2, 30.1, 30.0] + for _ in range(3) + ], + axis=0, +) +botm = [25.0, 10, 0.0] + +icelltype = 1 +K = 1000.0 +K33 = 100.0 +ss = 1e-6 +sy = 0.3 +strt_temp = 1.0 +scheme = "TVD" +dispersivity = 0.0 +ktw = 0.5918 +kts = 0.2700 +rhow = 1000.0 +cpw = 4183.0 +lhv = 2500.0 +cps = 760.0 +rhos = 1500.0 +K_therm_strmbed = [ + 1.5, + 1.75, + 2.0, +] # Thermal conductivity of the streambed material ($W/m/C$) +rbthcnd = 0.0001 +prsity = sy +drn_depth = 2.0 +ddrn = 1.0 +drn_cond = 1e5 + +ghbelv1 = 28.0 +ghbelv2 = 28.25 +ghbcond = 50000.0 +ghbtemp = 1.0 + +ghbspd = { + 0: [ + [0, 0, 0, ghbelv1, ghbcond, ghbtemp], + [0, 1, 0, ghbelv1, ghbcond, ghbtemp], + [0, 2, 0, ghbelv1, ghbcond, ghbtemp], + [1, 0, 0, ghbelv1, ghbcond, ghbtemp], + [1, 1, 0, ghbelv1, ghbcond, ghbtemp], + [1, 2, 0, ghbelv1, ghbcond, ghbtemp], + [2, 0, 0, ghbelv1, ghbcond, ghbtemp], + [2, 1, 0, ghbelv1, ghbcond, ghbtemp], + [2, 2, 0, ghbelv1, ghbcond, ghbtemp], + ], + 2: [ + [0, 0, 0, ghbelv2, ghbcond, ghbtemp], + [0, 1, 0, ghbelv2, ghbcond, ghbtemp], + [0, 2, 0, ghbelv2, ghbcond, ghbtemp], + [1, 0, 0, ghbelv2, ghbcond, ghbtemp], + [1, 1, 0, ghbelv2, ghbcond, ghbtemp], + [1, 2, 0, ghbelv2, ghbcond, ghbtemp], + [2, 0, 0, ghbelv2, ghbcond, ghbtemp], + [2, 1, 0, ghbelv2, ghbcond, ghbtemp], + [2, 2, 0, ghbelv2, ghbcond, ghbtemp], + ], +} + +# ifno cellid landflg ivertcn surfdp vks thtr thts thti eps [bndnm] +surfdep1 = 0.01 +surfdep2 = 0.001 +vks = 0.1 +eps = 4.0 +thti = 0.06 +thtr = 0.05 +thts = sy + thtr +iuzfbnd = np.array( + [ + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 1, 1, 1, 1, 0], + ] +) + +drn_pkdat = [] +uzf_pkdat = [] +uze_pkdat = [] +uze_perdat = [] +uzf_id_lkup = {} +ct = -1 +for k in np.arange(nlay): + landflg = 1 + if k > 0: + landflg = 0 + + for i in np.arange(nrow): + for j in np.arange(ncol): + if iuzfbnd[i, j] > 0: + ct += 1 + if k < nlay - 1: + ivertconn = ct + iuzfbnd.sum() + else: + ivertconn = -1 + + # gwf + uzf_pkdat.append( + [ + ct, + (k, i, j), + landflg, + ivertconn, + surfdep1, + vks, + thtr, + thts, + thti, + eps, + "uzf" + str(ct + 1).zfill(2), + ] + ) + # gwe + uze_pkdat.append((ct, 0)) + uze_perdat.append([ct, "INFILTRATION", 1.0]) + # generate a lookup dictionary based on the top layer + if k == 0: + drn_pkdat.append( + [(k, i, j), top[i, j] - drn_depth, drn_cond, ddrn] + ) + uzf_id_lkup.update({(i, j): ct}) + + +uze_perdat = { + 0: uze_perdat, +} + +mvr_pkdat = [] +for i in np.arange(nrow): + for j in np.arange(ncol - 1): + if iuzfbnd[i, j] > 0 and iuzfbnd[i, j + 1] > 0: + mvr_pkdat.append( + [ + "UZF-1", + uzf_id_lkup[(i, j)], + "UZF-1", + uzf_id_lkup[(i, j + 1)], + "FACTOR", + 1.0, + ] + ) + mvr_pkdat.append( + [ + "DRN-1", + uzf_id_lkup[(i, j)], + "UZF-1", + uzf_id_lkup[(i, j + 1)], + "FACTOR", + 1.0, + ] + ) + elif iuzfbnd[i, j] > 0 and j + 1 == ncol - 1: + mvr_pkdat.append( + ["UZF-1", uzf_id_lkup[(i, j)], "SFR-1", i, "FACTOR", 1.0] + ) + mvr_pkdat.append( + ["DRN-1", uzf_id_lkup[(i, j)], "SFR-1", i, "FACTOR", 1.0] + ) + +extdp = 3.0 +extwc = 0.05 +pet0 = 0.0 +pet1 = 0.001 +pet2 = 0.011 +finf0 = 0.2 +finf1 = 0.0 +finf2 = 0.2 +finf3 = 0.01 +zero = 0.0 +ntrail2 = 25 +nsets2 = 80 +uzf_spd = { + 0: [ + [0, finf0, pet0, extdp, extwc, zero, zero, zero], + [1, finf0, pet0, extdp, extwc, zero, zero, zero], + [2, finf0, pet0, extdp, extwc, zero, zero, zero], + [3, finf0, pet0, extdp, extwc, zero, zero, zero], + [4, finf0, pet0, extdp, extwc, zero, zero, zero], + [5, finf0, pet0, extdp, extwc, zero, zero, zero], + [6, finf0, pet0, extdp, extwc, zero, zero, zero], + [7, finf0, pet0, extdp, extwc, zero, zero, zero], + [8, finf1, pet0, extdp, extwc, zero, zero, zero], + [9, finf1, pet0, extdp, extwc, zero, zero, zero], + [10, finf1, pet0, extdp, extwc, zero, zero, zero], + [11, finf1, pet0, extdp, extwc, zero, zero, zero], + [12, finf1, pet0, extdp, extwc, zero, zero, zero], + [13, finf1, pet0, extdp, extwc, zero, zero, zero], + [14, finf1, pet0, extdp, extwc, zero, zero, zero], + [15, finf1, pet0, extdp, extwc, zero, zero, zero], + ], + 1: [], + 2: [], + 3: [ + [0, finf2, pet0, extdp, extwc, zero, zero, zero], + [1, finf2, pet0, extdp, extwc, zero, zero, zero], + [2, finf2, pet0, extdp, extwc, zero, zero, zero], + [3, finf2, pet0, extdp, extwc, zero, zero, zero], + [4, finf2, pet0, extdp, extwc, zero, zero, zero], + [5, finf2, pet0, extdp, extwc, zero, zero, zero], + [6, finf2, pet0, extdp, extwc, zero, zero, zero], + [7, finf2, pet0, extdp, extwc, zero, zero, zero], + [8, finf1, pet0, extdp, extwc, zero, zero, zero], + [9, finf1, pet0, extdp, extwc, zero, zero, zero], + [10, finf1, pet0, extdp, extwc, zero, zero, zero], + [11, finf1, pet0, extdp, extwc, zero, zero, zero], + [12, finf1, pet0, extdp, extwc, zero, zero, zero], + [13, finf1, pet0, extdp, extwc, zero, zero, zero], + [14, finf1, pet0, extdp, extwc, zero, zero, zero], + [15, finf1, pet0, extdp, extwc, zero, zero, zero], + ], +} + +rlen = delr +rwid = 1.0 +rgrd = 0.001 +rtp = 26.0 +rbth = 0.1 +rhk = 0.0 +rman = 0.02 +ncon = 2 +ustrf = 1.0 +ndv = 0 +pak_data = [] +for irno in range(nrow): + ncon = 2 + if irno in [0, nrow - 1]: + ncon = 1 + cellid = (0, irno, ncol - 1) + t = ( + irno, + cellid, + rlen, + rwid, + rgrd, + rtp, + rbth, + rhk, + rman, + ncon, + ustrf, + ndv, + ) + pak_data.append(t) + +con_data = [] +for irno in range(nrow): + if irno == 0: + t = (irno, -(irno + 1)) + elif irno == nrow - 1: + t = (irno, irno - 1) + else: + t = (irno, irno - 1, -(irno + 1)) + con_data.append(t) + +p_data = [ + (0, "INFLOW", 0.0), +] + + +def build_mf6_model(idx, ws): + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = cases[idx] + gwfname = "gwf-" + name + + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # create tdis package + flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, modelname=gwfname, newtonoptions="NEWTON", save_flows=True + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + complexity="MODERATE", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + sim.register_ims_package(ims, [gwf.name]) + + flopy.mf6.ModflowGwfdis( + gwf, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + ) + + # initial conditions + flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + flopy.mf6.ModflowGwfnpf(gwf, save_flows=True, icelltype=1, k=K, k33=K33) + + # aquifer storage + flopy.mf6.ModflowGwfsto(gwf, iconvert=1, ss=ss, sy=sy, transient=True) + + # ghb files + flopy.mf6.ModflowGwfghb( + gwf, + print_flows=True, + auxiliary="TEMPERATURE", + stress_period_data=ghbspd, + pname="GHB-1", + ) + + flopy.mf6.ModflowGwfdrn( + gwf, + save_flows=True, + print_input=True, + print_flows=True, + auxiliary="DDRN", + auxdepthname="DDRN", + stress_period_data=drn_pkdat, + mover=True, + pname="DRN-1", + ) + + # transient uzf info + uzf_obs = { + f"{gwfname}.uzfobs": [ + ("uzf01_dpth=0.5", "water-content", "uzf01", 0.5), + ( + "uzf01_dpth=1.5", + "water-content", + "uzf01", + 1.5, + ), # Relies on boundnames + ("uzf01_dpth=2.5", "water-content", "uzf01", 2.5), + ("uzf01_dpth=3.5", "water-content", "uzf01", 3.5), + ("uzf01_dpth=4.49", "water-content", "uzf01", 4.49), + ] + } + flopy.mf6.ModflowGwfuzf( + gwf, + mover=True, + print_flows=True, + save_flows=True, + wc_filerecord=gwfname + ".uzfwc.bin", + simulate_et=True, + simulate_gwseep=False, + linear_gwet=True, + boundnames=True, + observations=uzf_obs, + ntrailwaves=ntrail2, + nwavesets=nsets2, + nuzfcells=len(uzf_pkdat), + packagedata=uzf_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{gwfname}.uzf.bud", + pname="UZF-1", + filename=f"{gwfname}.uzf", + ) + + flopy.mf6.modflow.ModflowGwfsfr( + gwf, + save_flows=True, + print_input=True, + print_flows=True, + print_stage=True, + mover=True, + stage_filerecord=name + ".sfr.stg", + budget_filerecord=name + ".sfr.bud", + nreaches=nrow, + packagedata=pak_data, + pname="SFR-1", + connectiondata=con_data, + perioddata=p_data, + ) + + packages = [ + ("UZF-1",), + ("SFR-1",), + ("DRN-1",), + ] + flopy.mf6.ModflowGwfmvr( + gwf, + maxmvr=len(mvr_pkdat), + budget_filerecord=f"{gwfname}.mvr.bud", + maxpackages=len(packages), + print_flows=True, + packages=packages, + perioddata=mvr_pkdat, + ) + + # output control + 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")], + filename=f"{gwfname}.oc", + ) + + # ---------- + # GWE model + # ---------- + + gwename = "gwe-" + name + gwe = flopy.mf6.ModflowGwe( + sim, modelname=gwename, model_nam_file=f"{gwename}.nam" + ) + gwe.name_file.save_flows = True + + imsgwe = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + 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="{}.ims".format(gwename), + ) + sim.register_ims_package(imsgwe, [gwe.name]) + + # Instantiating MODFLOW 6 transport discretization package + flopy.mf6.ModflowGwedis( + gwe, + nogrb=True, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + pname="DIS", + filename=f"{gwename}.dis", + ) + + # Instantiating MODFLOW 6 transport initial concentrations + flopy.mf6.ModflowGweic( + gwe, + strt=strt_temp, + pname="IC", + filename=f"{gwename}.ic", + ) + + # Instantiating MODFLOW 6 transport advection package + flopy.mf6.ModflowGweadv( + gwe, scheme=scheme, pname="ADV", filename="{}.adv".format(gwename) + ) + + # Instantiating MODFLOW 6 transport dispersion package + flopy.mf6.ModflowGwecnd( + gwe, + xt3d_off=False, + alh=dispersivity, + ath1=dispersivity, + ktw=ktw * 86400, + kts=kts * 86400, + pname="CND", + filename=f"{gwename}.cnd", + ) + + # Instantiating MODFLOW 6 transport mass storage package + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + cps=cps, + rhos=rhos, + packagedata=[cpw, rhow, lhv], + pname="EST", + filename=f"{gwename}.est", + ) + + # Source-sink mixing package to support GHB input + srctype = "AUX" + auxname = "TEMPERATURE" + pname = ["GHB-1"] + # Inpput to SSM is: + sources = [[itm, srctype, auxname] for itm in pname] + + flopy.mf6.ModflowGwessm( + gwe, + sources=sources, + pname="SSM", + filename=f"{gwename}.ssm", + ) + + # Instantiate Streamflow Energy Transport package + sfepackagedata = [] + for irno in range(nrow): + t = (irno, strt_temp, K_therm_strmbed[irno], rbthcnd) + sfepackagedata.append(t) + + sfeperioddata = {0: [(0, "INFLOW", strt_temp)]} + + flopy.mf6.modflow.ModflowGwesfe( + gwe, + boundnames=False, + save_flows=True, + print_input=False, + print_flows=False, + print_temperature=True, + temperature_filerecord=gwename + ".sfe.bin", + budget_filerecord=gwename + ".sfe.bud", + packagedata=sfepackagedata, + reachperioddata=sfeperioddata, + flow_package_name="SFR-1", + pname="SFE-1", + filename="{}.sfe".format(gwename), + ) + + # Instantiating MODFLOW 6 unsaturated zone energy transport + flopy.mf6.ModflowGweuze( + gwe, + flow_package_name="UZF-1", + boundnames=False, + save_flows=True, + print_input=True, + print_flows=True, + print_temperature=True, + temperature_filerecord=gwename + ".uze.bin", + budget_filerecord=gwename + ".uze.bud", + packagedata=uze_pkdat, + uzeperioddata=uze_perdat, + pname="UZE-1", + filename=f"{gwename}.uze", + ) + + # mover transport package + fname = f"{gwename}.mve.bud" + mve = flopy.mf6.modflow.ModflowGwemve( + gwe, print_flows=True, budget_filerecord=fname + ) + + # Instantiate Output Control package for transport + flopy.mf6.ModflowGweoc( + gwe, + temperature_filerecord="{}.ucn".format(gwename), + budget_filerecord=f"{gwename}.bud", + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + temperatureprintrecord=[ + ("COLUMNS", 3, "WIDTH", 20, "DIGITS", 8, "GENERAL") + ], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename="{}.oc".format(gwename), + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename="{}.gwfgwe".format(gwename), + ) + + return sim + + +def build_models(idx, test): + # Start by building the MF6 model + sim = build_mf6_model(idx, test.workspace) + + return sim, None + + +def check_output(idx, test): + ws = test.workspace + name = cases[idx] + gwfname = "gwf-" + name + gwename = "gwe-" + name + + # Get the model budget items + fname = os.path.join(ws, gwfname + ".cbc") + assert os.path.isfile(fname) + modobj = flopy.utils.CellBudgetFile( + fname, precision="double", verbose=True + ) + + # Get the MVR results from GWF + fname = os.path.join(ws, gwfname + ".mvr.bud") + assert os.path.isfile(fname) + mvrobj = flopy.utils.CellBudgetFile( + fname, precision="double", verbose=True + ) + + # Get the MVE results from GWE + fname = os.path.join(ws, gwename + ".mve.bud") + assert os.path.isfile(fname) + mveobj = flopy.utils.CellBudgetFile( + fname, precision="double", verbose=False + ) + + ckstpkper = mveobj.get_kstpkper() + + drnbud = modobj.get_data(text="DRN") + drn2mvr = modobj.get_data(text="DRN-TO-MVR") + # uzf2mvr = modobj.get_data(text="UZF-GWD TO-MVR") + mvrdat = mvrobj.get_data(text="MOVER-FLOW") + mvedat = mveobj.get_data(text="MVE-FLOW") + + msg0 = "Accumulated cascading runoff is not as expected" + msg1 = ( + "Rejected infiltration being passed to MVR where it should not " + "be happening" + ) + msg2 = ( + "The accumulated cascading runoff that is finally passed to SFR " + "is not as expected" + ) + msg3 = "There should be no UZF -> SFR MVR flows in rows 2 and 3" + msg4 = ( + "There should be no cascading flow resulting from gw discharge " + "to land surface during the first stress period, but at least " + "one is not equal to 0.0" + ) + msg5 = ( + "GW discharge simulated with DRN from the right-most active DRN " + "cells should be 0 in the first stress period, but at least one " + "is not equal to 0.0" + ) + msg6 = ( + "GW discharge simulated by DRN for the right-most active DRN " + "cells should not be 0.0 in stress period 3, but at least one " + "is equal to 0.0" + ) + msg7 = "Accumulated cascaded energy is not accumulating as expected" + msg8 = ( + "Energy associated with rejected infiltration is being passed " + "to MVE where it should not be happening" + ) + msg9 = ( + "The accumulated cascading energy associated with runoff that " + "is finally passed to SFR is not as expected" + ) + msg10 = "There should be no UZE -> SFE MVE flows in rows 2 and 3" + msg11 = ( + "There should be no cascading energy flow resulting from gw " + "discharge to land surface, but at least one is not equal " + "to 0.0" + ) + msg12 = ( + "There should be no energy being passed by DRN from the " + "right-most active DRN cells during the first stress period, " + "but at least one is not equal to 0.0" + ) + msg13 = ( + "Energy accompanying gw discharge simulated by DRN for the " + "right-most active DRN cells should not be 0.0 in stress period " + "3, but at least one is equal to 0.0" + ) + + accum_runoff = 0.0 + accum_energy = 0.0 + for i, current_kstpkper in enumerate(ckstpkper): + mvr_rawdat = mvrobj.get_data(kstpkper=current_kstpkper) + mve_rawdat = mveobj.get_data(kstpkper=current_kstpkper) + if i == 0: + for idx, (itm, itm_e) in enumerate(zip(mvr_rawdat, mve_rawdat)): + # idx == 0 checks the UZF -> UZF MVR/MVE connections + if idx == 0: + for ct, val in enumerate(itm): + if ct < 7: + assert np.isclose( + itm[ct][-1], -1 * (finf0 - vks) + accum_runoff + ), msg0 + assert np.isclose( + itm_e[ct][-1], + (finf0 - vks) * rhow * cpw + accum_energy, + ), msg7 + accum_runoff += -1 * (finf0 - vks) + accum_energy += (finf0 - vks) * rhow * cpw + elif ct == 8: + accum_runoff += -1 * (finf0 - vks) + accum_energy += (finf0 - vks) * rhow * cpw + else: + assert itm[ct][-1] == 0, msg1 + assert itm_e[ct][-1] == 0, msg8 + + # idx == 1 checks the UZF -> SFR MVR connections + if idx == 1: + for ct, val in enumerate(itm): + if ct == 0: + assert np.isclose(itm[ct][-1], accum_runoff), msg2 + assert np.isclose( + itm_e[ct][-1], accum_energy + ), msg9 + else: + assert itm[ct][-1] == 0, msg3 + assert itm_e[ct][-1] == 0, msg10 + + # idx == 6 checks the DRN -> UZF MVR connections + # (groundwater discharge that cascades downhill) + if idx == 6: + for ct, val in enumerate(itm): + assert itm[ct][-1] == 0, msg4 + assert itm_e[ct][-1] == 0, msg11 + + # idx == 7 check the DRN -> SFR MVR connections + # (gw discharge from the lowest cell to SFR) + if idx == 7: + for ct, val in enumerate(itm): + assert itm[ct][-1] == 0, msg5 + assert itm_e[ct][-1] == 0, msg12 + + # Expect DRN -> MVR flows in the third stress period + if i == 2: + accum_runoff = 0 + accum_energy = 0 + for idx, (itm, itm_e) in enumerate(zip(mvr_rawdat, mve_rawdat)): + # idx == 0 checks the UZF -> UZF MVR connections + if idx == 0: + for ct, val in enumerate(itm): + if ct < 7: + assert np.isclose( + itm[ct][-1], -1 * (finf0 - vks) + accum_runoff + ), msg0 + assert np.isclose( + itm_e[ct][-1], + (finf0 - vks) * rhow * cpw + accum_energy, + ), msg7 + accum_runoff += -1 * (finf0 - vks) + accum_energy += (finf0 - vks) * rhow * cpw + elif ct == 8: + accum_runoff += -1 * (finf0 - vks) + accum_energy += (finf0 - vks) * rhow * cpw + else: + assert itm[ct][-1] == 0, msg1 + assert itm_e[ct][-1] == 0, msg8 + + # idx == 1 checks the UZF -> SFR MVR connections + if idx == 1: + for ct, val in enumerate(itm): + if ct == 0: + assert np.isclose(itm[ct][-1], accum_runoff), msg2 + assert np.isclose( + itm_e[ct][-1], accum_energy + ), msg9 + else: + assert itm[ct][-1] == 0, msg3 + assert itm_e[ct][-1] == 0, msg10 + + # idx == 6 checks the DRN -> UZF MVR connections + # (groundwater discharge that cascades downhill) + if idx == 6: + for ct, val in enumerate(itm): + assert itm[ct][-1] == 0, msg4 + assert itm_e[ct][-1] == 0, msg11 + + # idx == 7 check the DRN -> SFR MVR connections + # (gw discharge from the lowest cell to SFR) + if idx == 7: + for ct, val in enumerate(itm): + assert itm[ct][-1] < 0, msg6 + assert itm_e[ct][-1] > 0, msg13 + + print("Finished running checks") + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + targets=targets, + ) + test.run() diff --git a/autotest/test_gwe_uze00.py b/autotest/test_gwe_uze00.py index 9fb9b6336c1..a196d0514d5 100644 --- a/autotest/test_gwe_uze00.py +++ b/autotest/test_gwe_uze00.py @@ -338,8 +338,8 @@ def build_models(idx, test): ath1=dispersivity, ktw=0.5918 * 86400, kts=0.2700 * 86400, - pname="DSP", - filename=f"{gwename}.dsp", + pname="CND", + filename=f"{gwename}.cnd", ) # Instantiating MODFLOW 6 transport mass storage package @@ -353,8 +353,8 @@ def build_models(idx, test): cps=760.0, rhos=1500.0, packagedata=[cpw, rhow, lhv], - pname="MST", - filename=f"{gwename}.mst", + pname="EST", + filename=f"{gwename}.est", ) # Instantiating MODFLOW 6 constant temperature boundary condition at diff --git a/doc/mf6io/gwe/gwe.tex b/doc/mf6io/gwe/gwe.tex index 89eccf3b161..c908f5a85dd 100644 --- a/doc/mf6io/gwe/gwe.tex +++ b/doc/mf6io/gwe/gwe.tex @@ -145,6 +145,10 @@ \subsection{Unsaturated-Zone Energy Transport (UZE) Package} \subsection{Flow Model Interface (FMI) Package} \input{gwe/fmi} +\newpage +\subsection{Mover Energy Transport (MVE) Package} +\input{gwe/mve} + \newpage \subsection{Groundwater Energy Transport (GWE) Exchange} \input{gwe/gwe-gwe} diff --git a/doc/mf6io/gwe/mve.tex b/doc/mf6io/gwe/mve.tex new file mode 100644 index 00000000000..a931c3421dd --- /dev/null +++ b/doc/mf6io/gwe/mve.tex @@ -0,0 +1,18 @@ +Mover Energy Transport (MVE) Package information is read from the file that is specified by ``MVE6'' as the file type. Only one MVE Package can be specified for a GWE model. + +The MVE Package is used to route thermal energy according to flows from the GWF Water Mover (MVR) Package. The MVE Package must be activated by the user if the MVR Package is active in the corresponding the GWF Model. Flows from the GWF MVR Package must be available to the GWE model either through activation of a GWF-GWE Exchange or through specification of ``GWFMOVER'' in the PACKAGEDATA block of the GWE FMI Package. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-mve-options.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/gwe-mve-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwe-mve-example.dat} + diff --git a/doc/mf6io/mf6ivar/dfn/gwe-mve.dfn b/doc/mf6io/mf6ivar/dfn/gwe-mve.dfn new file mode 100644 index 00000000000..e67e76ba471 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/gwe-mve.dfn @@ -0,0 +1,106 @@ +# --------------------- gwe mve options --------------------- +# flopy subpackage mve_filerecord mve perioddata perioddata +# flopy parent_name_type parent_model_or_package MFModel/MFPackage + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'mover'} + +block options +name print_flows +type keyword +reader urword +optional true +longname print calculated flows to listing file +description REPLACE print_flows {'{#1}': 'lake'} + +block options +name save_flows +type keyword +reader urword +optional true +longname save lake flows to budget file +description REPLACE save_flows {'{#1}': 'lake'} + +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. + + diff --git a/doc/mf6io/mf6ivar/examples/gwe-mve-example.dat b/doc/mf6io/mf6ivar/examples/gwe-mve-example.dat new file mode 100644 index 00000000000..1b2835d13cf --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwe-mve-example.dat @@ -0,0 +1,6 @@ +BEGIN OPTIONS + PRINT_INPUT + PRINT_FLOWS + SAVE_FLOWS + BUDGET FILEOUT mygwemodel.mve.bud +END OPTIONS diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 4390cd2ff52..8821d4b6ad0 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -326,6 +326,14 @@ | 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 | MVE | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of mover information will be written to the listing file immediately after it is read. | +| GWE | MVE | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of lake 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 | MVE | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that lake flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | +| GWE | MVE | OPTIONS | BUDGET | KEYWORD | keyword to specify that record corresponds to the budget. | +| GWE | MVE | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. | +| GWE | MVE | OPTIONS | BUDGETFILE | STRING | name of the binary output file to write budget information. | +| GWE | MVE | OPTIONS | BUDGETCSV | KEYWORD | keyword to specify that record corresponds to the budget CSV. | +| GWE | MVE | 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 | 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. | diff --git a/doc/mf6io/mf6ivar/mf6ivar.py b/doc/mf6io/mf6ivar/mf6ivar.py index 4dba6f2ec62..aaf21a3c1d7 100644 --- a/doc/mf6io/mf6ivar/mf6ivar.py +++ b/doc/mf6io/mf6ivar/mf6ivar.py @@ -654,7 +654,6 @@ def write_appendix(blocks): f.write("\n\n\\hline\n\\end{longtable}\n\\label{table:blocks}\n\\normalsize\n") - def get_dfn_files(models): def is_sim_dfn(stem): return "sim" in stem diff --git a/doc/mf6io/mf6ivar/tex/appendixA.tex b/doc/mf6io/mf6ivar/tex/appendixA.tex index 7915b054671..273ed6a303c 100644 --- a/doc/mf6io/mf6ivar/tex/appendixA.tex +++ b/doc/mf6io/mf6ivar/tex/appendixA.tex @@ -100,6 +100,8 @@ GWE & LKE & PACKAGEDATA & yes \\ GWE & LKE & PERIOD & yes \\ \hline +GWE & MVE & OPTIONS & yes \\ +\hline GWE & MWE & OPTIONS & yes \\ GWE & MWE & PACKAGEDATA & yes \\ GWE & MWE & PERIOD & yes \\ diff --git a/doc/mf6io/mf6ivar/tex/gwe-mve-desc.tex b/doc/mf6io/mf6ivar/tex/gwe-mve-desc.tex new file mode 100644 index 00000000000..faae2d20ff4 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mve-desc.tex @@ -0,0 +1,23 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\item \texttt{PRINT\_INPUT}---keyword to indicate that the list of mover information will be written to the listing file immediately after it is read. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that the list of lake 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 lake flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. + +\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. + +\end{description} + diff --git a/doc/mf6io/mf6ivar/tex/gwe-mve-options.dat b/doc/mf6io/mf6ivar/tex/gwe-mve-options.dat new file mode 100644 index 00000000000..7b9ca2ed7b1 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-mve-options.dat @@ -0,0 +1,7 @@ +BEGIN OPTIONS + [PRINT_INPUT] + [PRINT_FLOWS] + [SAVE_FLOWS] + [BUDGET FILEOUT ] + [BUDGETCSV FILEOUT ] +END OPTIONS diff --git a/src/Exchange/exg-gwegwe.f90 b/src/Exchange/exg-gwegwe.f90 index b168d6c4b0c..6289e7e01ef 100644 --- a/src/Exchange/exg-gwegwe.f90 +++ b/src/Exchange/exg-gwegwe.f90 @@ -800,7 +800,7 @@ subroutine read_mvt(this, iout) ! for gwtmodel1 so that a call to save flows has an associated dis ! object. call mvt_cr(this%mvt, this%name, this%inmvt, iout, this%gwemodel1%fmi, & - this%gwemodel1%eqnsclfac, & + this%gwemodel1%eqnsclfac, this%gwemodel1%depvartype, & gwfmodelname1=this%gwfmodelname1, & gwfmodelname2=this%gwfmodelname2, & fmi2=this%gwemodel2%fmi) diff --git a/src/Exchange/exg-gwtgwt.f90 b/src/Exchange/exg-gwtgwt.f90 index c9e66e8d27c..e6559ce6f48 100644 --- a/src/Exchange/exg-gwtgwt.f90 +++ b/src/Exchange/exg-gwtgwt.f90 @@ -797,7 +797,7 @@ subroutine read_mvt(this, iout) ! for gwtmodel1 so that a call to save flows has an associated dis ! object. call mvt_cr(this%mvt, this%name, this%inmvt, iout, this%gwtmodel1%fmi, & - this%gwtmodel1%eqnsclfac, & + this%gwtmodel1%eqnsclfac, this%gwtmodel1%depvartype, & gwfmodelname1=this%gwfmodelname1, & gwfmodelname2=this%gwfmodelname2, & fmi2=this%gwtmodel2%fmi) diff --git a/src/Model/GroundWaterEnergy/gwe.f90 b/src/Model/GroundWaterEnergy/gwe.f90 index e387f77791c..8ed3e4394b1 100644 --- a/src/Model/GroundWaterEnergy/gwe.f90 +++ b/src/Model/GroundWaterEnergy/gwe.f90 @@ -70,7 +70,7 @@ module GweModule character(len=LENPACKAGETYPE), dimension(GWE_NBASEPKG) :: GWE_BASEPKG data GWE_BASEPKG/'DIS6 ', 'DISV6', 'DISU6', ' ', ' ', & ! 5 &'IC6 ', 'FMI6 ', 'EST6 ', 'ADV6 ', ' ', & ! 10 - &'CND6 ', 'SSM6 ', 'MVT6 ', 'OC6 ', ' ', & ! 15 + &'CND6 ', 'SSM6 ', 'MVE6 ', 'OC6 ', ' ', & ! 15 &'OBS6 ', ' ', ' ', ' ', ' ', & ! 20 &30*' '/ ! 50 diff --git a/src/Model/TransportModel/tsp-mvt.f90 b/src/Model/TransportModel/tsp-mvt.f90 index 17108e4f29a..c5b5a86f001 100644 --- a/src/Model/TransportModel/tsp-mvt.f90 +++ b/src/Model/TransportModel/tsp-mvt.f90 @@ -6,7 +6,8 @@ module TspMvtModule use KindModule, only: DP, I4B use ConstantsModule, only: LINELENGTH, MAXCHARLEN, DZERO, LENPAKLOC, & - DNODATA, LENPACKAGENAME, TABCENTER, LENMODELNAME + DNODATA, LENPACKAGENAME, TABCENTER, & + LENMODELNAME, LENVARNAME use SimModule, only: store_error use BaseDisModule, only: DisBaseType @@ -23,6 +24,7 @@ module TspMvtModule public :: mvt_cr type, extends(NumericalPackageType) :: TspMvtType + character(len=LENMODELNAME) :: gwfmodelname1 = '' !< name of model 1 character(len=LENMODELNAME) :: gwfmodelname2 = '' !< name of model 2 (set to modelname 1 for single model MVT) integer(I4B), pointer :: maxpackages !< max number of packages @@ -36,10 +38,13 @@ module TspMvtModule type(BudgetObjectType), pointer :: mvrbudobj => null() !< pointer to the water mover budget object character(len=LENPACKAGENAME), & dimension(:), pointer, contiguous :: paknames => null() !< array of package names + character(len=LENVARNAME) :: depvartype = '' ! ! -- table objects type(TableType), pointer :: outputtab => null() + contains + procedure :: mvt_df procedure :: mvt_ar procedure :: mvt_rp @@ -66,7 +71,7 @@ module TspMvtModule !> @brief Create a new mover transport object !< subroutine mvt_cr(mvt, name_model, inunit, iout, fmi1, eqnsclfac, & - gwfmodelname1, gwfmodelname2, fmi2) + depvartype, gwfmodelname1, gwfmodelname2, fmi2) ! -- dummy type(TspMvtType), pointer :: mvt character(len=*), intent(in) :: name_model @@ -74,15 +79,15 @@ subroutine mvt_cr(mvt, name_model, inunit, iout, fmi1, eqnsclfac, & integer(I4B), intent(in) :: iout type(TspFmiType), intent(in), target :: fmi1 real(DP), intent(in), pointer :: eqnsclfac !< governing equation scale factor + character(len=LENVARNAME), intent(in) :: depvartype !< dependent variable type ('concentration' or 'temperature') character(len=*), intent(in), optional :: gwfmodelname1 character(len=*), intent(in), optional :: gwfmodelname2 type(TspFmiType), intent(in), target, optional :: fmi2 -! ------------------------------------------------------------------------------ ! ! -- Create the object allocate (mvt) ! - ! -- create name and memory path + ! -- Create name and memory path call mvt%set_names(1, name_model, 'MVT', 'MVT') ! ! -- Allocate scalars @@ -95,12 +100,12 @@ subroutine mvt_cr(mvt, name_model, inunit, iout, fmi1, eqnsclfac, & mvt%fmi1 => fmi1 mvt%fmi2 => fmi1 ! - ! -- set pointers + ! -- Set pointers if (present(fmi2)) then mvt%fmi2 => fmi2 end if ! - ! -- set model names + ! -- Set model names if (present(gwfmodelname1)) then mvt%gwfmodelname1 = gwfmodelname1 end if @@ -108,12 +113,16 @@ subroutine mvt_cr(mvt, name_model, inunit, iout, fmi1, eqnsclfac, & mvt%gwfmodelname2 = gwfmodelname2 end if ! - ! -- create the budget object + ! -- Create the budget object call budgetobject_cr(mvt%budobj, 'TRANSPORT MOVER') ! ! -- Store pointer to governing equation scale factor mvt%eqnsclfac => eqnsclfac ! + ! -- Store pointer to labels associated with the current model so that the + ! package has access to the corresponding dependent variable type + mvt%depvartype = depvartype + ! ! -- Return return end subroutine mvt_cr @@ -121,27 +130,24 @@ end subroutine mvt_cr !> @brief Define mover transport object !< subroutine mvt_df(this, dis) - ! -- modules ! -- dummy class(TspMvtType) :: this class(DisBaseType), pointer, intent(in) :: dis - ! -- local ! -- formats character(len=*), parameter :: fmtmvt = & "(1x,/1x,'MVT -- MOVER TRANSPORT PACKAGE, VERSION 1, 4/15/2020', & &' INPUT READ FROM UNIT ', i0, //)" -! ------------------------------------------------------------------------------ ! - ! -- set pointer to dis + ! -- Set pointer to dis this%dis => dis ! - ! -- print a message identifying the MVT package. + ! -- Print a message identifying the MVT package. write (this%iout, fmtmvt) this%inunit ! ! -- Initialize block parser call this%parser%Initialize(this%inunit, this%iout) ! - ! -- initialize the budget table writer + ! -- Initialize the budget table writer call budget_cr(this%budget, this%memoryPath) ! ! -- Read mvt options @@ -156,7 +162,6 @@ end subroutine mvt_df !! Store a pointer to mvrbudobj, which contains the simulated water !! mover flows from either a gwf model MVR package or from a gwf-gwf !! exchange MVR package. - !! !< subroutine set_pointer_mvrbudobj(this, mvrbudobj) class(TspMvtType) :: this @@ -167,13 +172,10 @@ end subroutine set_pointer_mvrbudobj !> @brief Allocate and read mover-for-transport information !< subroutine mvt_ar(this) - ! -- modules ! -- dummy class(TspMvtType) :: this - ! -- locals -! ------------------------------------------------------------------------------ ! - ! -- setup the output table + ! -- Setup the output table call this%mvt_setup_outputtab() ! ! -- Return @@ -187,20 +189,17 @@ subroutine mvt_rp(this) use TdisModule, only: kper, kstp ! -- dummy class(TspMvtType) :: this - ! -- local - ! -- formats -! ------------------------------------------------------------------------------ ! ! -- At this point, the mvrbudobj is available to set up the mvt budobj if (kper * kstp == 1) then ! - ! -- if mvt is for a single model then point to fmi1 + ! -- If mvt is for a single model then point to fmi1 !cdl todo: this needs to be called from GwtGwtExg somehow for the 2 model case if (associated(this%fmi1, this%fmi2)) then call this%set_pointer_mvrbudobj(this%fmi1%mvrbudobj) end if ! - ! -- set up the mvt budobject + ! -- Set up the mvt budobject call this%mvt_scan_mvrbudobj() call this%mvt_setup_budobj() ! @@ -222,7 +221,6 @@ end subroutine mvt_rp !! hand side of the transport matrix equations. !< subroutine mvt_fc(this, cnew1, cnew2) - ! -- modules ! -- dummy class(TspMvtType) :: this real(DP), intent(in), dimension(:), contiguous, target :: cnew1 @@ -238,7 +236,6 @@ subroutine mvt_fc(this, cnew1, cnew2) real(DP), dimension(:), contiguous, pointer :: cnew type(TspFmiType), pointer :: fmi_pr !< pointer to provider model fmi package type(TspFmiType), pointer :: fmi_rc !< pointer to receiver model fmi package -! ------------------------------------------------------------------------------ ! ! -- Add mover QC terms to the receiver packages nbudterm = this%mvrbudobj%nbudterm @@ -249,7 +246,8 @@ subroutine mvt_fc(this, cnew1, cnew2) ! -- Set pointers to the fmi packages for the provider and the receiver call this%set_fmi_pr_rc(i, fmi_pr, fmi_rc) ! - ! -- Set a pointer to the GWT model concentration associated with the provider + ! -- Set a pointer to the GWT model concentration (or temperature) + ! associated with the provider cnew => cnew1 if (associated(fmi_pr, this%fmi2)) then cnew => cnew2 @@ -280,8 +278,9 @@ subroutine mvt_fc(this, cnew1, cnew2) cp = DZERO if (fmi_pr%iatp(ipr) /= 0) then ! - ! -- Provider package is being represented by an APT (SFT, LKT, MWT, UZT) - ! so set the concentration to the simulated concentration of APT + ! -- Provider package is being represented by an APT (SFT, LKT, MWT, UZT, + ! SFE, LKE, MWE, UZE) so set the dependent variable (concentration or + ! temperature) to the simulated dependent variable of the APT. cp = concpak(id1) else ! @@ -290,10 +289,10 @@ subroutine mvt_fc(this, cnew1, cnew2) ! SFT, LKT, MWT, or UZT, so use the GWT cell concentration igwtnode = fmi_pr%gwfpackages(ipr)%nodelist(id1) cp = cnew(igwtnode) - + ! end if ! - ! -- add the mover rate times the provider concentration into the receiver + ! -- Add the mover rate times the provider concentration into the receiver ! make sure these are accumulated since multiple providers can move ! water into the same receiver if (fmi_rc%iatp(irc) /= 0) then @@ -323,40 +322,40 @@ subroutine set_fmi_pr_rc(this, ibudterm, fmi_pr, fmi_rc) integer(I4B), intent(in) :: ibudterm type(TspFmiType), pointer :: fmi_pr type(TspFmiType), pointer :: fmi_rc - + ! fmi_pr => null() fmi_rc => null() if (this%gwfmodelname1 == '' .and. this%gwfmodelname2 == '') then fmi_pr => this%fmi1 fmi_rc => this%fmi1 else - ! modelname for provider is this%mvrbudobj%budterm(i)%text1id1 + ! -- Modelname for provider is this%mvrbudobj%budterm(i)%text1id1 if (this%mvrbudobj%budterm(ibudterm)%text1id1 == this%gwfmodelname1) then - ! -- model 1 is the provider + ! -- Model 1 is the provider fmi_pr => this%fmi1 else if (this%mvrbudobj%budterm(ibudterm)%text1id1 == & this%gwfmodelname2) then - ! -- model 2 is the provider + ! -- Model 2 is the provider fmi_pr => this%fmi2 else - ! must be an error + ! -- Must be an error !cdl todo: programming error print *, this%mvrbudobj%budterm(ibudterm)%text1id1 print *, this%gwfmodelname1 print *, this%gwfmodelname2 stop "error in set_fmi_pr_rc" end if - - ! modelname for receiver is this%mvrbudobj%budterm(i)%text1id2 + ! + ! -- Modelname for receiver is this%mvrbudobj%budterm(i)%text1id2 if (this%mvrbudobj%budterm(ibudterm)%text1id2 == this%gwfmodelname1) then - ! -- model 1 is the receiver + ! -- Model 1 is the receiver fmi_rc => this%fmi1 else if (this%mvrbudobj%budterm(ibudterm)%text1id2 == & this%gwfmodelname2) then - ! -- model 2 is the receiver + ! -- Model 2 is the receiver fmi_rc => this%fmi2 else - ! must be an error + ! -- Must be an error !cdl todo: programming error print *, this%mvrbudobj%budterm(ibudterm)%text1id2 print *, this%gwfmodelname1 @@ -364,7 +363,7 @@ subroutine set_fmi_pr_rc(this, ibudterm, fmi_pr, fmi_rc) stop "error in set_fmi_pr_rc" end if end if - + ! if (.not. associated(fmi_pr)) then print *, 'Could not find FMI Package...' stop "error in set_fmi_pr_rc" @@ -388,7 +387,6 @@ subroutine mvt_cc(this, kiter, iend, icnvgmod, cpak, dpak) integer(I4B), intent(in) :: icnvgmod character(len=LENPAKLOC), intent(inout) :: cpak real(DP), intent(inout) :: dpak - ! -- local ! -- formats character(len=*), parameter :: fmtmvrcnvg = & "(/,1x,'MOVER PACKAGE REQUIRES AT LEAST TWO OUTER ITERATIONS. CONVERGE & @@ -403,25 +401,22 @@ subroutine mvt_cc(this, kiter, iend, icnvgmod, cpak, dpak) end if end if ! - ! -- return + ! -- Return return end subroutine mvt_cc !> @brief Write mover terms to listing file !< subroutine mvt_bd(this, cnew1, cnew2) - ! -- modules ! -- dummy class(TspMvtType) :: this real(DP), dimension(:), contiguous, intent(in) :: cnew1 real(DP), dimension(:), contiguous, intent(in) :: cnew2 - ! -- local -! ------------------------------------------------------------------------------ ! - ! -- fill the budget object + ! -- Fill the budget object call this%mvt_fill_budobj(cnew1, cnew2) ! - ! -- return + ! -- Return return end subroutine mvt_bd @@ -436,7 +431,6 @@ subroutine mvt_ot_saveflow(this, icbcfl, ibudfl) integer(I4B), intent(in) :: ibudfl ! -- locals integer(I4B) :: ibinun -! ------------------------------------------------------------------------------ ! ! -- Save the mover flows from the budobj to a mover binary file ibinun = 0 @@ -456,13 +450,10 @@ end subroutine mvt_ot_saveflow !> @brief Print mover flow table !< subroutine mvt_ot_printflow(this, icbcfl, ibudfl) - ! -- modules ! -- dummy class(TspMvtType) :: this integer(I4B), intent(in) :: icbcfl integer(I4B), intent(in) :: ibudfl - ! -- locals -! ------------------------------------------------------------------------------ ! ! -- Print the mover flow table if (ibudfl /= 0 .and. this%iprflow /= 0) then @@ -494,28 +485,21 @@ subroutine mvt_ot_bdsummary(this, ibudfl) ! ! -- Accumulate the rates do i = 1, this%maxpackages - do j = 1, this%budobj%nbudterm - do n = 1, this%budobj%budterm(j)%nlist - ! - ! -- provider is inflow to mover + ! -- Provider is inflow to mover if (this%paknames(i) == this%budobj%budterm(j)%text2id1) then ratin(i) = ratin(i) + this%budobj%budterm(j)%flow(n) end if ! - ! -- receiver is outflow from mover + ! -- Receiver is outflow from mover if (this%paknames(i) == this%budobj%budterm(j)%text2id2) then ratout(i) = ratout(i) + this%budobj%budterm(j)%flow(n) end if - end do - end do - end do - ! ! -- Send rates to budget object call this%budget%reset() @@ -554,25 +538,23 @@ subroutine mvt_da(this) use MemoryManagerModule, only: mem_deallocate ! -- dummy class(TspMvtType) :: this - ! -- local -! ------------------------------------------------------------------------------ ! ! -- Deallocate arrays if package was active if (this%inunit > 0) then ! - ! -- character array + ! -- Character array deallocate (this%paknames) ! - ! -- budget object + ! -- Budget object call this%budget%budget_da() deallocate (this%budget) ! - ! -- budobj + ! -- Budobj call this%budobj%budgetobject_da() deallocate (this%budobj) nullify (this%budobj) ! - ! -- output table object + ! -- Output table object if (associated(this%outputtab)) then call this%outputtab%table_da() deallocate (this%outputtab) @@ -588,7 +570,7 @@ subroutine mvt_da(this) call mem_deallocate(this%ibudgetout) call mem_deallocate(this%ibudcsv) ! - ! -- deallocate scalars in NumericalPackageType + ! -- Deallocate scalars in NumericalPackageType call this%NumericalPackageType%da() ! ! -- Return @@ -604,10 +586,8 @@ subroutine allocate_scalars(this) use MemoryManagerModule, only: mem_allocate, mem_setptr ! -- dummy class(TspMvtType) :: this - ! -- local -! ------------------------------------------------------------------------------ ! - ! -- allocate scalars in NumericalPackageType + ! -- Allocate scalars in NumericalPackageType call this%NumericalPackageType%allocate_scalars() ! ! -- Allocate @@ -637,18 +617,18 @@ subroutine read_options(this) character(len=MAXCHARLEN) :: fname integer(I4B) :: ierr logical :: isfound, endOfBlock + ! -- formats character(len=*), parameter :: fmtflow = & "(4x, a, 1x, a, 1x, ' WILL BE SAVED TO FILE: ', a, & &/4x, 'OPENED ON UNIT: ', I0)" character(len=*), parameter :: fmtflow2 = & &"(4x, 'FLOWS WILL BE SAVED TO BUDGET FILE')" -! ------------------------------------------------------------------------------ ! - ! -- get options block + ! -- Get options block call this%parser%GetBlock('OPTIONS', isfound, ierr, blockRequired=.false., & supportOpenClose=.true.) ! - ! -- parse options block if detected + ! -- Parse options block if detected if (isfound) then write (this%iout, '(1x,a)') 'PROCESSING MVT OPTIONS' do @@ -702,7 +682,7 @@ subroutine read_options(this) write (this%iout, '(1x,a)') 'END OF MVT OPTIONS' end if ! - ! -- return + ! -- Return return end subroutine read_options @@ -722,15 +702,18 @@ subroutine mvt_setup_budobj(this) character(len=LENMODELNAME) :: modelname1, modelname2 character(len=LENPACKAGENAME) :: packagename1, packagename2 character(len=LENBUDTXT) :: text -! ------------------------------------------------------------------------------ ! ! -- Assign terms to set up the mover budget object nbudterm = this%mvrbudobj%nbudterm ncv = 0 - text = ' MVT-FLOW' naux = 0 + if (this%depvartype == 'CONCENTRATION') then + text = ' MVT-FLOW' + else + text = ' MVE-FLOW' + end if ! - ! -- set up budobj + ! -- Set up budobj call this%budobj%budgetobject_df(ncv, nbudterm, 0, 0, bddim_opt='M') ! ! -- Go through the water mover budget terms and set up the transport @@ -757,7 +740,6 @@ end subroutine mvt_setup_budobj !> @brief Copy mover-for-transport flow terms into this%budobj !< subroutine mvt_fill_budobj(this, cnew1, cnew2) - ! -- modules ! -- dummy class(TspMvtType) :: this real(DP), intent(in), dimension(:), contiguous, target :: cnew1 @@ -777,8 +759,6 @@ subroutine mvt_fill_budobj(this, cnew1, cnew2) real(DP) :: cp real(DP) :: q real(DP) :: rate - ! -- formats -! ----------------------------------------------------------------------------- ! ! -- Go through the water mover budget terms and set up the transport ! mover budget terms @@ -813,7 +793,7 @@ subroutine mvt_fill_budobj(this, cnew1, cnew2) rate = -q * cp * this%eqnsclfac end if ! - ! -- add the rate to the budterm + ! -- Add the rate to the budterm call this%budobj%budterm(i)%update_term(n1, n2, rate) end do end do @@ -821,7 +801,7 @@ subroutine mvt_fill_budobj(this, cnew1, cnew2) ! --Terms are filled, now accumulate them for this time step call this%budobj%accumulate_terms() ! - ! -- return + ! -- Return return end subroutine mvt_fill_budobj @@ -848,13 +828,13 @@ subroutine mvt_scan_mvrbudobj(this) end do this%maxpackages = maxpackages ! - ! -- allocate paknames + ! -- Allocate paknames allocate (this%paknames(this%maxpackages)) do i = 1, this%maxpackages this%paknames(i) = '' end do ! - ! -- scan through mvrbudobj and create unique paknames + ! -- Scan through mvrbudobj and create unique paknames ipos = 1 do i = 1, nbudterm found = .false. @@ -886,14 +866,14 @@ subroutine mvt_setup_outputtab(this) integer(I4B) :: maxrow integer(I4B) :: ilen ! - ! -- allocate and initialize the output table + ! -- Allocate and initialize the output table if (this%iprflow /= 0) then ! - ! -- dimension table + ! -- Dimension table ntabcol = 7 maxrow = 0 ! - ! -- initialize the output table object + ! -- Initialize the output table object title = 'TRANSPORT MOVER PACKAGE ('//trim(this%packName)// & ') FLOW RATES' call table_cr(this%outputtab, this%packName, title) @@ -915,10 +895,10 @@ subroutine mvt_setup_outputtab(this) call this%outputtab%initialize_column(text, ilen) text = 'RECEIVER ID' call this%outputtab%initialize_column(text, 10) - + ! end if ! - ! -- return + ! -- Return return end subroutine mvt_setup_outputtab @@ -937,16 +917,15 @@ subroutine mvt_print_outputtab(this) integer(I4B) :: inum integer(I4B) :: ntabrows integer(I4B) :: nlist -! ------------------------------------------------------------------------------ ! - ! -- determine number of table rows + ! -- Determine number of table rows ntabrows = 0 do i = 1, this%budobj%nbudterm nlist = this%budobj%budterm(i)%nlist ntabrows = ntabrows + nlist end do ! - ! -- set table kstp and kper + ! -- Set table kstp and kper call this%outputtab%set_kstpkper(kstp, kper) ! ! -- Add terms and print the table @@ -975,7 +954,7 @@ subroutine mvt_print_outputtab(this) end do end do ! - ! -- return + ! -- Return return end subroutine mvt_print_outputtab diff --git a/src/Model/TransportModel/tsp.f90 b/src/Model/TransportModel/tsp.f90 index e40f3abc0ba..18a2bc20c1e 100644 --- a/src/Model/TransportModel/tsp.f90 +++ b/src/Model/TransportModel/tsp.f90 @@ -791,7 +791,7 @@ subroutine create_tsp_packages(this, indis) mempathic = mempath case ('FMI6') this%infmi = inunit - case ('MVT6') + case ('MVT6', 'MVE6') this%inmvt = inunit case ('ADV6') this%inadv = inunit @@ -816,7 +816,7 @@ subroutine create_tsp_packages(this, indis) call ssm_cr(this%ssm, this%name, this%inssm, this%iout, this%fmi, & this%eqnsclfac, this%depvartype) call mvt_cr(this%mvt, this%name, this%inmvt, this%iout, this%fmi, & - this%eqnsclfac) + this%eqnsclfac, this%depvartype) call oc_cr(this%oc, this%name, this%inoc, this%iout) call tsp_obs_cr(this%obs, this%inobs) ! diff --git a/src/Utilities/BudgetObject.f90 b/src/Utilities/BudgetObject.f90 index 0e60fd605a6..c0f935fa9ce 100644 --- a/src/Utilities/BudgetObject.f90 +++ b/src/Utilities/BudgetObject.f90 @@ -123,7 +123,7 @@ subroutine budgetobject_df(this, ncv, nbudterm, iflowja, nsto, & character(len=16) :: labeltitle character(len=20) :: bdzone ! - ! -- set values + ! -- Set values this%ncv = ncv this%nbudterm = nbudterm this%iflowja = iflowja @@ -288,7 +288,7 @@ subroutine accumulate_terms(this) use TdisModule, only: delt ! -- dummy class(BudgetObjectType) :: this - ! -- dummy + ! -- local character(len=LENBUDTXT) :: flowtype integer(I4B) :: i real(DP) :: ratin, ratout @@ -526,7 +526,6 @@ subroutine read_flows(this, dis, ibinun) real(DP) :: delt real(DP) :: pertim real(DP) :: totim - ! -- dummy integer(I4B) :: i ! ! -- Read flows for each budget term @@ -544,7 +543,7 @@ end subroutine read_flows subroutine budgetobject_da(this) ! -- dummy class(BudgetObjectType) :: this - ! -- dummy + ! -- local integer(I4B) :: i ! ! -- Save flows for each budget term @@ -658,8 +657,9 @@ subroutine bfr_advance(this, dis, iout) class(BudgetObjectType) :: this class(DisBaseType), intent(in) :: dis integer(I4B), intent(in) :: iout - ! -- dummy + ! -- local logical :: readnext + ! -- formats character(len=*), parameter :: fmtkstpkper = & &"(1x,/1x, a, ' READING BUDGET TERMS FOR KSTP ', i0, ' KPER ', i0)" character(len=*), parameter :: fmtbudkstpkper = & @@ -708,7 +708,7 @@ subroutine fill_from_bfr(this, dis, iout) class(BudgetObjectType) :: this class(DisBaseType), intent(in) :: dis integer(I4B), intent(in) :: iout - ! -- dummy + ! -- local integer(I4B) :: i logical :: success !