diff --git a/autotest/test_gwe_uze00.py b/autotest/test_gwe_uze00.py new file mode 100644 index 00000000000..5ef6a0e885e --- /dev/null +++ b/autotest/test_gwe_uze00.py @@ -0,0 +1,612 @@ +# +# - Outer columns not active for unsaturated zone, but are present to host +# constant head boundaries at the bottom of the model. +# +# +-------+ +# |///////| = Inactive cell +# +-------+ +# +# Model depiction: +# +# +-------+-------+-------+ +# |///////| |///////| Layer 1 +# +-------+-------+-------+ +# |///////| |///////| Layer 2 +# +-------+-------+-------+ +# |///////| |///////| Layer 3 +# +-------+-------+-------+ +# |///////| |///////| +# + -- -- + -- -- + -- -- + +# |///////| |///////| Layer x (Middle portion of model not shown) +# + -- -- + -- -- + -- -- + +# |///////| |///////| +# +-------+-------+-------+ +# | | | | Layer 99 +# +-------+-------+-------+ +# | | | | Layer 100 +# +-------+-------+-------+ + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework + +import flopy.utils.binaryfile as bf +import math + +import matplotlib.pyplot as plt + + +# Analytical solution, from Barends (2010) Equation 5 +def temp_analyt(t, z, t0, tinfil, v, d): + if t == 0.0: + temp = t0 + else: + denom = 2.0 * math.sqrt(d * t) + ztermm = (z - v * t) / denom + ztermp = (z + v * t) / denom + vterm = v * z / d + if vterm < 100.0: + # might need to adjust this limit + temp = t0 + 0.5 * (tinfil - t0) * ( + math.erfc(ztermm) + math.exp(vterm) * math.erfc(ztermp) + ) + else: + zeta = 1.0 / (1.0 + 0.47047 * ztermp) + polyterm = zeta * ( + 0.3480242 + zeta * (-0.0958798 + zeta * 0.7478556) + ) + temp = t0 + 0.5 * (tinfil - t0) * ( + math.erfc(ztermm) + math.exp(vterm - ztermp ** 2) * polyterm + ) + + return temp + + +# Model units +length_units = "meters" +time_units = "days" + +nlay, nrow, ncol = 101, 1, 3 +nper = 2 +perlen = [1.0e9, 100.0] +nstp = [1, 100] +tsmult = len(perlen) * [1.0] + +delr = 1.0 +delc = 1.0 +delz = 0.1 # 10 cm +strt = 0.05 +top = 10.0005 +botm = [ + 9.9995 +] # Top layer is very thin for application of the boundary condition +for i in np.arange(1, nlay): + bot = 10.0 - (i * delz) + botm.append(round(bot, 1)) + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-9, 1e-3, 0.97 +steady = {0: False, 1: False} +transient = {0: True, 1: True} + +idomain_u = [0, 1, 0] +idomain_l = [1, 1, 1] +idomain = [] +for i in np.arange(nlay): + if i < 99: + idomain.append(idomain_u) + else: + idomain.append(idomain_l) + +idomain = np.array(idomain) + +strt_temp = 10.0 +scheme = "UPSTREAM" +dispersivity = 0.0 +prsity = 0.2 + +# transient uzf info +# iuzno cellid landflg ivertcn surfdp vks thtr thts thti eps [bndnm] +uzf_pkdat = [[0, (0, 0, 1), 1, 1, 0.00001, 1, 0.0001, 0.20, 0.055, 4]] + +# Continue building the UZF list of objects +for iuzno in np.arange(1, 101, 1): + if iuzno < 99: + ivertconn = iuzno + 1 + else: + ivertconn = -1 + + uzf_pkdat.append( + [iuzno, (iuzno, 0, 1), 0, ivertconn, 0.01, 1, 0.0001, 0.20, 0.055, 4] + ) + +iuz_cell_dict = {} +cell_iuz_dict = {} +for i, itm in enumerate(uzf_pkdat): + iuz_cell_dict.update({itm[0]: (itm[1][0], itm[1][1], itm[1][2])}) + cell_iuz_dict.update({(itm[1][0], itm[1][1], itm[1][2]): itm[0]}) + +finf = 0.01 +extdp = 0.0 +pet = 0.0 +extwc = 0.0 +zero = 0.0 +uzf_spd = { + 0: [[0, finf, pet, extdp, extwc, zero, zero, zero]], + 1: [[0, finf, pet, extdp, extwc, zero, zero, zero]], +} + +cases = ["uze00"] + + +def build_models(idx, test): + name = cases[idx] + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # create tdis package + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + flopy.mf6.ModflowTdis( + sim, time_units=time_units, nper=nper, perioddata=tdis_rc + ) + + gwfname = "gwf_" + name + gwename = "gwe_" + name + + newtonoptions = ["NEWTON", "UNDER_RELAXATION"] + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + newtonoptions=newtonoptions, + 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, + filename=f"{gwfname}.ims", + ) + sim.register_ims_package(ims, [gwf.name]) + + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=idomain, + filename=f"{gwfname}.dis", + ) + + # initial conditions + flopy.mf6.ModflowGwfic( + gwf, + strt=strt, + filename=f"{gwfname}.ic", + ) + + # node property flow + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=100.0, + k33=10, + filename=f"{gwfname}.npf", + ) + + # aquifer storage + flopy.mf6.ModflowGwfsto( + gwf, + iconvert=1, + ss=1e-5, + sy=prsity, + steady_state=steady, + transient=transient, + filename=f"{gwfname}.sto", + ) + + # chd files + chdval = 0.05 + chdspd = {0: [[(100, 0, 0), chdval, 10.0], [(100, 0, 2), chdval, 10.0]]} + + flopy.mf6.ModflowGwfchd( + gwf, + auxiliary=["TEMPERATURE"], + print_flows=True, + stress_period_data=chdspd, + pname="CHD-1", + filename=f"{gwfname}.chd", + ) + + # Unsaturated-zone flow package + flopy.mf6.ModflowGwfuzf( + gwf, + print_flows=True, + save_flows=True, + wc_filerecord=gwfname + ".uzfwc.bin", + simulate_et=False, + simulate_gwseep=False, + linear_gwet=False, + boundnames=False, + ntrailwaves=15, + nwavesets=40, + nuzfcells=len(uzf_pkdat), + packagedata=uzf_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{gwfname}.uzf.bud", + pname="UZF-1", + filename=f"{gwfname}.uzf", + ) + + # output control + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwfname}.oc", + ) + + # ---------------------------------- + # Instantiating MODFLOW 6 GWE model + # ---------------------------------- + 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, + idomain=idomain, + 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=0.5918 * 86400, + kts=0.2700 * 86400, + pname="DSP", + filename=f"{gwename}.dsp", + ) + + # Instantiating MODFLOW 6 transport mass storage package + rhow = 1000.0 + cpw = 4183.0 + lhv = 2500.0 + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + cps=760.0, + rhos=1500.0, + packagedata=[cpw, rhow, lhv], + pname="MST", + filename=f"{gwename}.mst", + ) + + # Instantiating MODFLOW 6 constant temperature boundary condition at + ctpspd = {0: [[(0, 0, 1), 10.0]], 1: [[(0, 0, 1), 20.0]]} + + flopy.mf6.ModflowGwectp( + gwe, + save_flows=True, + print_flows=True, + stress_period_data=ctpspd, + pname="CTP", + filename="{}.ctp".format(gwename), + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + srctype = "AUX" + auxname = "TEMPERATURE" + pname = ["CHD-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", + ) + + # Instantiating MODFLOW 6 energy transport source-sink mixing package + uzepackagedata = [(iuz, 10.0) for iuz in range(nlay)] + uzeperioddata = { + 0: [[0, "INFILTRATION", 10.0]], + 1: [[0, "INFILTRATION", 20.0]], + } + + 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=uzepackagedata, + uzeperioddata=uzeperioddata, + pname="UZE-1", + filename=f"{gwename}.uze", + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGweoc( + gwe, + pname="OC", + budget_filerecord="{}.cbc".format(gwename), + temperature_filerecord="{}.ucn".format(gwename), + temperatureprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename="{}.gwfgwe".format(gwename), + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating flow...") + + name = cases[idx] + gwfname = "gwf_" + name + gwename = "gwe_" + name + ws = test.workspace + + # check some output... + wc_fl = gwfname + ".uzfwc.bin" + wcobj = flopy.utils.HeadFile(os.path.join(ws, wc_fl), text="water-content") + wc = wcobj.get_alldata() + + fl2 = gwename + ".uze.bin" + + uzeobj = flopy.utils.HeadFile(os.path.join(ws, fl2), text="TEMPERATURE") + temps = uzeobj.get_alldata() + + t = np.linspace(0.0, 100.0, 101) + z = np.arange(0.05, 10.0, 0.1) + z = np.insert(z, 0, 0.0) + + t0 = 10.0 + tinfil = 20.0 + + q = finf # infiltration rate + rhos = 1500.0 + Cps = 760.0 + rhow = 1000.0 + Cpw = 4183.0 + rhowCpw = Cpw * rhow + rhosCps = Cps * rhos + + Kts = 23328.0 + Ktw = 0.0 + + steady_wc = wc[1, 0, 0, 1] + Sw = steady_wc / prsity + + rhoCp_bulk = Sw * prsity * rhowCpw + (1 - prsity) * rhosCps + Kt_bulk = Sw * prsity * Ktw + (1 - prsity) * Kts + v = rhowCpw / rhoCp_bulk * q + D = Kt_bulk / rhoCp_bulk + + # Put analytical solution in place + analytical_sln = np.zeros((len(t), len(z))) + for i, tm in enumerate(t): + for j, depth in enumerate(z): + temp = temp_analyt(tm, depth, t0, tinfil, v, D) + analytical_sln[i, j] = temp + + # Run checks + msg0 = ( + "Simulated solution no longer falling within" + " default tolerance where it previously did" + ) + # Compare day 1. For layer 20 and below, the defaults of allclose should work + assert np.allclose(analytical_sln[1, 19:], temps[1, 0, 0, 19:]), msg0 + # Compare day 10. For layer 39 and below, the defaults of allclose should work + assert np.allclose(analytical_sln[10, 38:], temps[10, 0, 0, 38:]), msg0 + # Compare day 50. For layer 84 and below, the defaults of allclose should work + assert np.allclose(analytical_sln[50, 83:], temps[50, 0, 0, 83:]), msg0 + # Compare day 100, fits are generally good, but do not pass allclose default settings + + # Ensure that the differences in the 1st day fall within established bounds + msg1 = ( + "Simulated fits to analytical solution are " + "falling outside established bounds on day 1" + ) + assert ( + np.max(analytical_sln[1, :18] - temps[1, 0, 0, :18]) <= 1.52921097880 + ), msg1 + assert ( + np.min(analytical_sln[1, :18] - temps[1, 0, 0, :18]) >= -0.32260871278 + ), msg1 + + # Ensure that the differences on day 10 fall within established bounds + msg2 = ( + "Simulated fits to analytical solution are " + "falling outside established bounds on day 10" + ) + assert ( + np.max(analytical_sln[10, :37] - temps[10, 0, 0, :37]) <= 0.15993441016 + ), msg2 + assert ( + np.min(analytical_sln[10, :37] - temps[10, 0, 0, :37]) + >= -0.22298707253 + ), msg2 + + # Ensure that the differences on day 50 fall within established bounds + msg3 = ( + "Simulated fits to analytical solution are " + "falling outside established bounds on day 50" + ) + assert ( + np.max(analytical_sln[50, :82] - temps[50, 0, 0, :82]) <= 0.09327747258 + ), msg3 + assert ( + np.min(analytical_sln[50, :82] - temps[50, 0, 0, :82]) + >= -0.21182907402 + ), msg3 + + # Ensure that the differences on day 50 fall within established bounds + msg3 = ( + "Simulated fits to analytical solution are " + "falling outside established bounds on day 50" + ) + assert ( + np.max(analytical_sln[50, :82] - temps[50, 0, 0, :82]) <= 0.09327747258 + ), msg3 + assert ( + np.min(analytical_sln[50, :82] - temps[50, 0, 0, :82]) + >= -0.21182907402 + ), msg3 + + # Ensure that the differences on day 100 fall within established bounds + msg4 = ( + "Simulated fits to analytical solution are " + "falling outside established bounds on day 100" + ) + assert np.max(analytical_sln[100] - temps[100]) <= 0.10680304268, msg4 + assert np.min(analytical_sln[100] - temps[100]) >= -0.20763221276, msg4 + + # If a plot is needed for visual inspection, change following if statement to "True" + if False: + analytical_sln = np.zeros((len(t), len(z))) + for i, tm in enumerate(t): + for j, depth in enumerate(z): + temp = temp_analyt(tm, depth, t0, tinfil, v, D) + analytical_sln[i, j] = temp + + # first transient stress period + line1 = plt.plot( + analytical_sln[1], z, "-", color="red", label="Analytical" + ) + line2 = plt.plot( + temps[1, 0, 0], z, "-.", color="blue", label="MODFLOW 6" + ) + # 10th transient stress period + plt.plot(analytical_sln[10], z, "-", color="red") + plt.plot(temps[10, 0, 0], z, "-.", color="blue") + # 50th transient stress period + plt.plot(analytical_sln[50], z, "-", color="red") + plt.plot(temps[50, 0, 0], z, "-.", color="blue") + # last stress period + plt.plot(analytical_sln[100], z, "-", color="red") + plt.plot(temps[100, 0, 0], z, "-.", color="blue") + # add labels + plt.text(11.0, 0.85, "1 day", fontsize=10) + plt.text(12.0, 1.65, "10 days", fontsize=10) + plt.text(14.0, 2.90, "50 days", fontsize=10) + plt.text(16.0, 4.00, "100 days", fontsize=10) + + plt.gca().invert_yaxis() + plt.xlabel( + "$Temperature, C$" + ) # For latex replace with: '$Temperature, ^{\circ}C$' + plt.ylabel("$Depth, m$") + plt.minorticks_on() + plt.axhline(y=0.0) + plt.legend(loc="lower right", frameon=False) + plt.savefig(os.path.join(ws, "fit_view.png"), format="png") + + +# - 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/autotest/test_gwe_uze00_flux.py b/autotest/test_gwe_uze00_flux.py new file mode 100644 index 00000000000..3873865d2a6 --- /dev/null +++ b/autotest/test_gwe_uze00_flux.py @@ -0,0 +1,753 @@ +# - Similar to test_gwe_uze00.py; however, this looks into whether the +# flux input is correct without pinning the temperature of the top-most +# cell to a specified value. +# +# - Outer columns not active for unsaturated zone, but are present to host +# constant head boundaries at the bottom of the model. +# +# +-------+ +# |///////| = Inactive cell +# +-------+ +# +# Model depiction: +# +# +-------+-------+-------+ +# |///////| |///////| Layer 1 +# +-------+-------+-------+ +# |///////| |///////| Layer 2 +# +-------+-------+-------+ +# |///////| |///////| Layer 3 +# +-------+-------+-------+ +# |///////| |///////| +# + -- -- + -- -- + -- -- + +# |///////| |///////| Layer x (Middle portion of model not shown) +# + -- -- + -- -- + -- -- + +# |///////| |///////| +# +-------+-------+-------+ +# | | | | Layer 99 +# +-------+-------+-------+ +# | | | | Layer 100 +# +-------+-------+-------+ + +import os + +import flopy +import numpy as np +import pytest +import math + +import flopy.utils.binaryfile as bf +from framework import TestFramework + + +# Analytical solution derived by Alden, similar in form to +# Barends (2010) Equation 5 - but remember that that solution +# pins the temperature of the top cell to a specified temperature +def flux_analyt(t, z, qt0, qtinfil, v, d): + if t == 0.0: + flux = qt0 + else: + denom = 2.0 * math.sqrt(d * t) + ztermm = (z - v * t) / denom + ztermp = (z + v * t) / denom + vterm = v * z / d + if vterm < 100.0: + # might need to adjust this limit + flux = qt0 + (qtinfil - qt0) * 0.5 * ( + math.erfc(ztermm) + math.exp(vterm) * math.erfc(ztermp) + ) + else: + zeta = 1.0 / (1.0 + 0.47047 * ztermp) + polyterm = zeta * ( + 0.3480242 + zeta * (-0.0958798 + zeta * 0.7478556) + ) + flux = qt0 + 0.5 * (qtinfil - qt0) * ( + math.erfc(ztermm) + math.exp(vterm - ztermp ** 2) * polyterm + ) + return flux + + +def temp_analyt(t, z, t0, tinfil, v, d): + if t == 0.0: + temp = t0 + else: + denom = 2.0 * math.sqrt(d * t) + ztermm = (z - v * t) / denom + ztermp = (z + v * t) / denom + vterm = v * z / d + if vterm < 100.0: + # might need to adjust this limit + temp = t0 + 0.5 * (tinfil - t0) * ( + math.erfc(ztermm) + math.exp(vterm) * math.erfc(ztermp) + ) + else: + zeta = 1.0 / (1.0 + 0.47047 * ztermp) + polyterm = zeta * ( + 0.3480242 + zeta * (-0.0958798 + zeta * 0.7478556) + ) + temp = t0 + 0.5 * (tinfil - t0) * ( + math.erfc(ztermm) + math.exp(vterm - ztermp ** 2) * polyterm + ) + return temp + + +# Model units +length_units = "meters" +time_units = "days" + +nlay, nrow, ncol = 101, 1, 3 +nper = 2 +perlen = [1.0e9, 100.0] +nstp = [1, 100] +tsmult = len(perlen) * [1.0] + +delr = 1.0 +delc = 1.0 +delz = 0.1 # 10 cm +strt = 0.05 +top = 10.0005 +botm = [ + 9.9995 +] # Top layer is very thin for application of the boundary condition +for i in np.arange(1, nlay): + bot = 10.0 - (i * delz) + botm.append(round(bot, 1)) + +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-9, 1e-3, 0.97 +steady = {0: False, 1: False} +transient = {0: True, 1: True} + +idomain_u = [0, 1, 0] +idomain_l = [1, 1, 1] +idomain = [] +for i in np.arange(nlay): + if i < 99: + idomain.append(idomain_u) + else: + idomain.append(idomain_l) + +idomain = np.array(idomain) + +strt_temp = 10.0 +scheme = "UPSTREAM" +dispersivity = 0.0 +prsity = 0.2 + +# transient uzf info +# iuzno cellid landflg ivertcn surfdp vks thtr thts thti eps [bndnm] +uzf_pkdat = [[0, (0, 0, 1), 1, 1, 0.00001, 1, 0.0001, 0.20, 0.055, 4]] + +# Continue building the UZF list of objects +for iuzno in np.arange(1, 101, 1): + if iuzno < nlay - 1: + ivertconn = iuzno + 1 + else: + ivertconn = -1 + + uzf_pkdat.append( + [iuzno, (iuzno, 0, 1), 0, ivertconn, 0.01, 1, 0.0001, 0.20, 0.055, 4] + ) + +iuz_cell_dict = {} +cell_iuz_dict = {} +for i, itm in enumerate(uzf_pkdat): + iuz_cell_dict.update({itm[0]: (itm[1][0], itm[1][1], itm[1][2])}) + cell_iuz_dict.update({(itm[1][0], itm[1][1], itm[1][2]): itm[0]}) + +finf = 0.01 +extdp = 0.0 +pet = 0.0 +extwc = 0.0 +zero = 0.0 +uzf_spd = { + 0: [[0, finf, pet, extdp, extwc, zero, zero, zero]], + 1: [[0, finf, pet, extdp, extwc, zero, zero, zero]], +} + +cases = ["uze00_flux"] + + +def build_models(idx, test): + name = cases[idx] + + # build MODFLOW 6 files + ws = test.workspace + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + + # create tdis package + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + flopy.mf6.ModflowTdis( + sim, time_units=time_units, nper=nper, perioddata=tdis_rc + ) + + gwfname = "gwf_" + name + gwename = "gwe_" + name + + newtonoptions = ["NEWTON", "UNDER_RELAXATION"] + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=gwfname, + newtonoptions=newtonoptions, + 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, + filename=f"{gwfname}.ims", + ) + sim.register_ims_package(ims, [gwf.name]) + + flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=idomain, + filename=f"{gwfname}.dis", + ) + + # initial conditions + flopy.mf6.ModflowGwfic( + gwf, + strt=strt, + filename=f"{gwfname}.ic", + ) + + # node property flow + flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=100.0, + k33=10, + filename=f"{gwfname}.npf", + ) + + # aquifer storage + flopy.mf6.ModflowGwfsto( + gwf, + iconvert=1, + ss=1e-5, + sy=prsity, + steady_state=steady, + transient=transient, + filename=f"{gwfname}.sto", + ) + + # chd files + chdval = 0.05 + chdspd = {0: [[(100, 0, 0), chdval, 10.0], [(100, 0, 2), chdval, 10.0]]} + + flopy.mf6.ModflowGwfchd( + gwf, + auxiliary=["TEMPERATURE"], + print_flows=True, + stress_period_data=chdspd, + pname="CHD-1", + filename=f"{gwfname}.chd", + ) + + # Unsaturated-zone flow package + flopy.mf6.ModflowGwfuzf( + gwf, + print_flows=True, + save_flows=True, + wc_filerecord=gwfname + ".uzfwc.bin", + simulate_et=False, + simulate_gwseep=False, + linear_gwet=False, + boundnames=False, + ntrailwaves=15, + nwavesets=40, + nuzfcells=len(uzf_pkdat), + packagedata=uzf_pkdat, + perioddata=uzf_spd, + budget_filerecord=f"{gwfname}.uzf.bud", + pname="UZF-1", + filename=f"{gwfname}.uzf", + ) + + # output control + flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwfname}.oc", + ) + + # ---------------------------------- + # Instantiating MODFLOW 6 GWE model + # ---------------------------------- + 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, + idomain=idomain, + 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=0.5918 * 86400, + kts=0.2700 * 86400, + pname="DSP", + filename=f"{gwename}.dsp", + ) + + # Instantiating MODFLOW 6 transport mass storage package + rhow = 1000.0 + cpw = 4183.0 + lhv = 2500.0 + flopy.mf6.ModflowGweest( + gwe, + save_flows=True, + porosity=prsity, + cps=760.0, + rhos=1500.0, + packagedata=[cpw, rhow, lhv], + pname="MST", + filename=f"{gwename}.mst", + ) + + # Instantiating MODFLOW 6 transport source-sink mixing package + srctype = "AUX" + auxname = "TEMPERATURE" + pname = ["CHD-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", + ) + + # Instantiating MODFLOW 6 energy transport source-sink mixing package + uzepackagedata = [(iuz, 10.0) for iuz in range(nlay)] + uzeperioddata = { + 0: [[0, "INFILTRATION", 10.0]], + 1: [[0, "INFILTRATION", 20.0]], + } + + 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=uzepackagedata, + uzeperioddata=uzeperioddata, + pname="UZE-1", + filename=f"{gwename}.uze", + ) + + # Instantiate MODFLOW 6 heat transport output control package + flopy.mf6.ModflowGweoc( + gwe, + pname="OC", + budget_filerecord="{}.cbc".format(gwename), + temperature_filerecord="{}.ucn".format(gwename), + temperatureprintrecord=[ + ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + printrecord=[("TEMPERATURE", "ALL"), ("BUDGET", "ALL")], + filename=f"{gwename}.oc", + ) + + # Instantiate Gwf-Gwe Exchange package + flopy.mf6.ModflowGwfgwe( + sim, + exgtype="GWF6-GWE6", + exgmnamea=gwfname, + exgmnameb=gwename, + filename="{}.gwfgwe".format(gwename), + ) + + return sim, None + + +def check_output(idx, test): + print("evaluating flow...") + + name = cases[idx] + gwfname = "gwf_" + name + gwename = "gwe_" + name + ws = test.workspace + + # check some output... + wc_fl = gwfname + ".uzfwc.bin" + wcobj = flopy.utils.HeadFile(os.path.join(ws, wc_fl), text="water-content") + wc = wcobj.get_alldata() + + # temperature output + fl2 = gwename + ".uze.bin" + uzeobj = flopy.utils.HeadFile(os.path.join(ws, fl2), text="TEMPERATURE") + temps = uzeobj.get_alldata() + + # Cell flows output + qfile = gwename + ".cbc" + gweflowsobj = flopy.utils.CellBudgetFile(os.path.join(ws, qfile)) + + # Binary grid file needed for post-processing + fgrb = gwfname + ".dis.grb" + grb_file = os.path.join(ws, fgrb) + + # UZE flows + fuzebud = gwename + ".uze.bud" + uzeflowsobj = flopy.utils.CellBudgetFile(os.path.join(ws, fuzebud)) + flowsadv = uzeflowsobj.get_data(text="FLOW-JA-FACE") + + t = np.linspace(0.0, 100.0, 101) + z = np.linspace(0.0, 9.9, 99) + + q = finf # infiltration rate + area = delr * delc + rhos = 1500.0 + Cps = 760.0 + rhow = 1000.0 + Cpw = 4183.0 + rhowCpw = Cpw * rhow + rhosCps = Cps * rhos + + Kts = 23328.0 + Ktw = 0.0 + + steady_wc = wc[1, 0, 0, 1] + Sw = steady_wc / prsity + + rhoCp_bulk = Sw * prsity * rhowCpw + (1 - prsity) * rhosCps + Kt_bulk = Sw * prsity * Ktw + (1 - prsity) * Kts + v = rhowCpw / rhoCp_bulk * q + D = Kt_bulk / rhoCp_bulk + + t0 = 10.0 + tinfil = 20.0 + qt0 = q * t0 + qtinfil = q * tinfil + + # for converting from J/day to W + unitadj = 1 / 86400 + + # Get analytical solution + conv10 = [] + cond10 = [] + conv50 = [] + cond50 = [] + conv100 = [] + cond100 = [] + analytical_sln = np.zeros((len(t), len(z))) + simulated_sln = np.zeros((len(t), len(z))) + for i, tm in enumerate(t): + if i == 0: + gweflowjaface = gweflowsobj.get_data( + text="FLOW-JA-FACE", kstpkper=(0, 0) + ) + else: + gweflowjaface = gweflowsobj.get_data( + text="FLOW-JA-FACE", kstpkper=(i - 1, 1) + ) + flowscond = flopy.mf6.utils.postprocessing.get_structured_faceflows( + gweflowjaface[0][0], grb_file=grb_file + ) + for j, depth in enumerate(z): + fluxa = flux_analyt(tm, depth, qt0, qtinfil, v, D) + analytical_sln[i, j] = fluxa * rhowCpw * unitadj + + (uze1, uze2, floadv) = flowsadv[i][2 * j + 1] + (fjunk1, flocond, fjunk2) = flowscond[1][j][0] + flo = floadv * rhowCpw * unitadj + flocond * rhowCpw * unitadj + if i == 10: + conv10.append(floadv * rhowCpw * unitadj) + cond10.append(flocond * rhowCpw * unitadj) + elif i == 50: + conv50.append(floadv * rhowCpw * unitadj) + cond50.append(flocond * rhowCpw * unitadj) + elif i == 100: + conv100.append(floadv * rhowCpw * unitadj) + cond100.append(flocond * rhowCpw * unitadj) + + flux = flo / area + simulated_sln[i, j] = flux + + # Run checks + msg0 = ( + "Simulated solution has deviated too far from the analytical solution" + ) + # Following values are calculated as "percent differences" when determining if fits are acceptable + # Day 10 + assert ( + np.max( + ((analytical_sln[10] * rhowCpw) - simulated_sln[10]) + / (analytical_sln[10] * rhowCpw) + * 100 + ) + <= 0.51091366512 + ), msg0 + assert ( + np.min( + ((analytical_sln[10] * rhowCpw) - simulated_sln[10]) + / (analytical_sln[10] * rhowCpw) + * 100 + ) + >= -2.62119104308 + ), msg0 + + # Day 50 + assert ( + np.max( + ((analytical_sln[50] * rhowCpw) - simulated_sln[50]) + / (analytical_sln[50] * rhowCpw) + * 100 + ) + <= 0.3171034702 + ), msg0 + assert ( + np.min( + ((analytical_sln[50] * rhowCpw) - simulated_sln[50]) + / (analytical_sln[50] * rhowCpw) + * 100 + ) + >= -2.52020856408 + ), msg0 + + # Day 100 + assert ( + np.max( + ((analytical_sln[100] * rhowCpw) - simulated_sln[100]) + / (analytical_sln[100] * rhowCpw) + * 100 + ) + <= 0.37655443171 + ), msg0 + assert ( + np.min( + ((analytical_sln[100] * rhowCpw) - simulated_sln[100]) + / (analytical_sln[100] * rhowCpw) + * 100 + ) + >= -2.56004275226 + ), msg0 + + # If plot is needed, change next statement to "if True:" + if False: + import matplotlib.pyplot as plt + from matplotlib import transforms + from matplotlib.collections import PathCollection + from matplotlib.patches import Patch + from matplotlib.lines import Line2D + + fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(10, 5)) + # 10 days + # ------- + polys1 = ax1.stackplot( + z, + conv10, + cond10, + labels=["Convection", "Conduction"], + colors=["lightseagreen", "lightgreen"], + ) + ax1.set_xlim((-0.05, 10.05)) + xlims = ax1.get_xlim() + for poly in polys1: + for path in poly.get_paths(): + path.vertices = path.vertices[:, ::-1] + ax1.set_xlim(0.095 * rhowCpw * unitadj, 0.205 * rhowCpw * unitadj) + ax1.set_ylim(xlims[::-1]) + ax1.plot(simulated_sln[10], z, "-", color="blue", linewidth=1) + ax1.plot(analytical_sln[10], z, "-.", color="red") + ax1.text(4.9, 0.4, "10 Days") + + legend_elements = [ + Line2D( + [0], + [0], + linestyle="-", + color="blue", + lw=1, + label="MODFLOW 6 Total Heat Flux", + ), + Line2D( + [0], [0], linestyle="-.", color="red", lw=1, label="Analytical" + ), + Patch( + facecolor="lightseagreen", + edgecolor="lightseagreen", + label="Convection", + ), + Patch( + facecolor="lightgreen", + edgecolor="lightgreen", + label="Conduction", + ), + ] + + ax1.legend(handles=legend_elements, loc="lower right", frameon=False) + ax1.set_xlabel("Energy Flux, $\dfrac{Watts}{m^2}$") + ax1.set_ylabel("Depth, m") + + # 50 days + # ------- + polys2 = ax2.stackplot( + z, + conv50, + cond50, + labels=["Convection", "Conduction"], + colors=["lightseagreen", "lightgreen"], + ) + ax2.set_xlim((-0.05, 10.05)) + xlims = ax2.get_xlim() + for poly in polys2: + for path in poly.get_paths(): + path.vertices = path.vertices[:, ::-1] + ax2.set_xlim(0.095 * rhowCpw * unitadj, 0.205 * rhowCpw * unitadj) + ax2.set_ylim(xlims[::-1]) + ax2.plot(simulated_sln[50], z, "-", color="blue", linewidth=1) + ax2.plot(analytical_sln[50], z, "-.", color="red") + ax2.set_xlabel("Energy Flux, $\dfrac{Watts}{m^2}$") + ax2.text(4.9, 0.4, "50 Days") + + # 100 days + # ------- + polys3 = ax3.stackplot( + z, + conv100, + cond100, + labels=["Convection", "Conduction"], + colors=["lightseagreen", "lightgreen"], + ) + ax3.set_xlim((-0.05, 10.05)) + xlims = ax3.get_xlim() + for poly in polys3: + for path in poly.get_paths(): + path.vertices = path.vertices[:, ::-1] + ax3.set_xlim(0.095 * rhowCpw * unitadj, 0.205 * rhowCpw * unitadj) + ax3.set_ylim(xlims[::-1]) + ax3.plot(simulated_sln[100], z, "-", color="blue", linewidth=1) + ax3.plot(analytical_sln[100], z, "-.", color="red") + ax3.set_xlabel("Energy Flux, $\dfrac{Watts}{m^2}$") + ax3.text(4.9, 0.4, "100 Days") + + plt.tight_layout() + plt.savefig(os.path.join(ws, "dual_view.png"), format="png") + + line1 = plt.plot( + analytical_sln[10], z, "-", color="red", label="Analytical" + ) + line2 = plt.plot( + simulated_sln[10], z, "-.", color="blue", label="MODFLOW 6" + ) + # 50th transient stress period + plt.plot(analytical_sln[50], z, "-", color="red") + plt.plot(simulated_sln[50], z, "-.", color="blue") + # last stress period + plt.plot(analytical_sln[100], z, "-", color="red") + plt.plot(simulated_sln[100], z, "-.", color="blue") + # add labels + plt.text(11.0, 0.85, "1 day", fontsize=10) + plt.text(12.0, 1.65, "10 days", fontsize=10) + plt.text(14.0, 2.90, "50 days", fontsize=10) + plt.text(16.0, 4.00, "100 days", fontsize=10) + + plt.gca().invert_yaxis() + plt.xlabel("$Energy Flux, -$") + plt.ylabel("$Depth, m$") + plt.minorticks_on() + plt.axhline(y=0.0) + plt.legend(loc="lower right", frameon=False) + plt.savefig(os.path.join(ws, "fit_view.png"), format="png") + + +# - 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-uzeobs.tex b/doc/Common/gwe-uzeobs.tex index f9a6aa1f2ba..241c86ef653 100644 --- a/doc/Common/gwe-uzeobs.tex +++ b/doc/Common/gwe-uzeobs.tex @@ -1,14 +1,14 @@ % general APT observations -UZE & temperature & ifno or boundname & -- & uze cell temperature. If boundname is specified, boundname must be unique for each uze cell. \\ -UZE & flow-ja-face & ifno or boundname & ifno or -- & Energy flow between two uze cells. If a boundname is specified for ID1, then the result is the total energy flow for all uze cells. If a boundname is specified for ID1 then ID2 is not used.\\ -UZE & storage & ifno or boundname & -- & Simulated energy storage flow rate for a uze cell or group of uze cells. \\ -UZE & constant & ifno or boundname & -- & Simulated energy constant-flow rate for a uze cell or a group of uze cells. \\ -UZE & from-mvr & ifno or boundname & -- & Simulated energy inflow into a uze cell or group of uze cells from the MVE package. Energy inflow is calculated as the product of provider temperature and the mover flow rate. \\ -UZE & uze & ifno or boundname & -- & Energy flow rate for a uze cell or group of uze cells and its aquifer connection(s). \\ +UZE & temperature & uzeno or boundname & -- & uze cell temperature. If boundname is specified, boundname must be unique for each uze cell. \\ +UZE & flow-ja-face & uzeno or boundname & uzeno or -- & Energy flow between two uze cells. If a boundname is specified for ID1, then the result is the total energy flow for all uze cells. If a boundname is specified for ID1 then ID2 is not used.\\ +UZE & storage & uzeno or boundname & -- & Simulated energy storage flow rate for a uze cell or group of uze cells. \\ +UZE & constant & uzeno or boundname & -- & Simulated energy constant-flow rate for a uze cell or a group of uze cells. \\ +UZE & from-mvr & uzeno or boundname & -- & Simulated energy inflow into a uze cell or group of uze cells from the MVE package. Energy inflow is calculated as the product of provider temperature and the mover flow rate. \\ +UZE & uze & uzeno or boundname & -- & Energy flow rate for a uze cell or group of uze cells and its aquifer connection(s). \\ %observations specific to the uze package % infiltration rej-inf uzet rej-inf-to-mvr -UZE & infiltration & ifno or boundname & -- & Infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature. \\ -UZE & rej-inf & ifno or boundname & -- & Rejected infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature. \\ -UZE & uzet & ifno or boundname & -- & Unsaturated zone evapotranspiration rate applied to a uze cell or group of uze cells multiplied by the latent heat of vaporization, heat capacity of the simulated fluid, and density of the simulated flued. \\ -UZE & rej-inf-to-mvr & ifno or boundname & -- & Rejected infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature that is sent to the mover package. \\ +UZE & infiltration & uzeno or boundname & -- & Infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature. \\ +UZE & rej-inf & uzeno or boundname & -- & Rejected infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature. \\ +UZE & uzet & uzeno or boundname & -- & Unsaturated zone evapotranspiration rate applied to a uze cell or group of uze cells multiplied by the uze cell temperature. \\ +UZE & rej-inf-to-mvr & uzeno or boundname & -- & Rejected infiltration rate applied to a uze cell or group of uze cells multiplied by the infiltration temperature that is sent to the mover package. \\ diff --git a/doc/mf6io/gwe/gwe.tex b/doc/mf6io/gwe/gwe.tex index 0e7e30211bc..c4ab64fe833 100644 --- a/doc/mf6io/gwe/gwe.tex +++ b/doc/mf6io/gwe/gwe.tex @@ -137,6 +137,10 @@ \subsection{Lake Energy Transport (LKE) Package} \subsection{Multi-Aquifer Well Energy Transport (MWE) Package} \input{gwe/mwe} +\newpage +\subsection{Unsaturated-Zone Energy Transport (UZE) Package} +\input{gwe/uze} + \newpage \subsection{Flow Model Interface (FMI) Package} \input{gwe/fmi} diff --git a/doc/mf6io/gwe/namefile.tex b/doc/mf6io/gwe/namefile.tex index 4419af99530..d8cc1d18fec 100644 --- a/doc/mf6io/gwe/namefile.tex +++ b/doc/mf6io/gwe/namefile.tex @@ -37,6 +37,7 @@ \subsubsection{Explanation of Variables} SFE6 & Streamflow Energy Transport Package & * \\ LKE6 & Lake Energy Transport Package & * \\ MWE6 & Multi-Aquifer Well Energy Transport Package & * \\ +UZE6 & Unsaturated-Zone Energy Transport Package & * \\ OBS6 & Observations Option \\ \hline \end{tabular*} diff --git a/doc/mf6io/gwe/uze.tex b/doc/mf6io/gwe/uze.tex new file mode 100644 index 00000000000..1e65d5bbd11 --- /dev/null +++ b/doc/mf6io/gwe/uze.tex @@ -0,0 +1,55 @@ +Unsaturated Zone Energy Transport (UZE) Package information is read from the file that is specified by ``UZE6'' as the file type. There can be as many UZE Packages as necessary for a GWE model. Each UZE Package is designed to work with flows from a corresponding GWF UZF Package. By default \mf uses the UZE package name to determine which UZF Package corresponds to the UZE Package. Therefore, the package name of the UZE Package (as specified in the GWE name file) must match with the name of the corresponding UZF 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 UZE Package cannot be used without a corresponding GWF UZF Package. + +The UZE Package does not have a dimensions block; instead, dimensions for the UZE Package are set using the dimensions from the corresponding UZF Package. For example, the UZF Package requires specification of the number of cells (NUZFCELLS). UZE sets the number of UZE cells equal to NUZFCELLS. Therefore, the PACKAGEDATA block below must have NUZFCELLS entries in it. + +\vspace{5mm} +\subsubsection{Structure of Blocks} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-uze-options.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-uze-packagedata.dat} +\lstinputlisting[style=blockdefinition]{./mf6ivar/tex/gwe-uze-period.dat} + +\vspace{5mm} +\subsubsection{Explanation of Variables} +\begin{description} +\input{./mf6ivar/tex/gwe-uze-desc.tex} +\end{description} + +\vspace{5mm} +\subsubsection{Example Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwe-uze-example.dat} + +\vspace{5mm} +\subsubsection{Available observation types} +Unsaturated Zone Energy Transport Package observations include UZF cell temperature and all of the terms that contribute to the continuity equation for each UZF cell. Additional UZE Package observations include energy flow rates for individual UZF cells, or groups of UZF cells. The data required for each UZE Package observation type is defined in table~\ref{table:gwe-uzeobstype}. Negative and positive values for \texttt{uzt} 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 UZE package, respectively. + +\begin{longtable}{p{2cm} p{2.75cm} p{2cm} p{1.25cm} p{7cm}} +\caption{Available UZE 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 UZE 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-uzeobs.tex} +\label{table:gwe-uzeobstype} +\end{longtable} + +\vspace{5mm} +\subsubsection{Example Observation Input File} +\lstinputlisting[style=inputfile]{./mf6ivar/examples/gwe-uze-example-obs.dat} + + diff --git a/doc/mf6io/mf6ivar/dfn/gwe-uze.dfn b/doc/mf6io/mf6ivar/dfn/gwe-uze.dfn new file mode 100644 index 00000000000..1f272617b73 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/gwe-uze.dfn @@ -0,0 +1,438 @@ +# --------------------- gwe uze 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 concentration 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 concentrations 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}': 'unsaturated zone flow'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'unsaturated zone flow'} + +block options +name print_temperature +type keyword +reader urword +optional true +longname print calculated temperatures to listing file +description REPLACE print_temperature {'{#1}': 'UZF cell', '{#2}': 'temperatures', '{#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}': 'unsaturated zone'} + +block options +name save_flows +type keyword +reader urword +optional true +longname save UZE cell flows to budget file +description REPLACE save_flows {'{#1}': 'unsaturated zone'} + +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}': 'UZE'} + + +# --------------------- gwe uze packagedata --------------------- + +block packagedata +name packagedata +type recarray uzfno strt aux boundname +shape (maxbound) +reader urword +longname +description + +block packagedata +name uzfno +type integer +shape +tagged false +in_record true +reader urword +longname UZF cell number for this entry +description integer value that defines the UZF cell number associated with the specified PACKAGEDATA data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. Unsaturated zone flow information must be specified for every UZF cell or the program will terminate with an error. The program also will terminate with an error if information for a UZF cell 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 UZF cell temperature +description real value that defines the starting temperature for the unsaturated zone flow cell. + +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}': 'unsaturated zone flow'} + +block packagedata +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname UZF cell name +description REPLACE boundname {'{#1}': 'unsaturated zone flow'} + + +# --------------------- gwe uze 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 uzeperioddata +type recarray uzfno uzesetting +shape +reader urword +longname +description + +block period +name uzfno +type integer +shape +tagged false +in_record true +reader urword +longname unsaturated zone flow cell number for this entry +description integer value that defines the UZF cell number associated with the specified PERIOD data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. +numeric_index true + +block period +name uzesetting +type keystring status temperature infiltration uzet 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 UZESETTING string include: STATUS, TEMPERATURE, INFILTRATION, UZET, and AUXILIARY. These settings are used to assign the temperature associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. + +block period +name status +type string +shape +tagged true +in_record true +reader urword +longname unsaturated zone flow cell temperature status +description keyword option to define UZF cell status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the UZF cell. If a UZF cell is inactive, then there will be no energy fluxes into or out of the UZF cell and the inactive value will be written for the UZF cell temperature. If a UZF cell is constant, then the temperature for the UZF cell 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 unsaturated zone flow cell temperature +description real or character value that defines the temperature for the unsaturated zone flow cell. The specified TEMPERATURE is only applied if the unsaturated zone flow cell is a constant temperature cell. 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 infiltration +type string +shape +tagged true +in_record true +reader urword +time_series true +longname infiltration temperature +description real or character value that defines the temperature of the infiltration $(^\circ C)$ for the UZF cell. 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 uzet +type string +shape +tagged true +in_record true +reader urword +time_series true +longname unsaturated zone et temperature +description real or character value that states what fraction of the simulated unsaturated zone evapotranspiration is associated with evaporation. The evaporative losses from the unsaturated zone moisture content will have an evaporative cooling effect on the GWE cell from which the evaporation originated. If this value is larger than 1, then it will be reset to 1. 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-uze-example-obs.dat b/doc/mf6io/mf6ivar/examples/gwe-uze-example-obs.dat new file mode 100644 index 00000000000..3859fd40ff1 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwe-uze-example-obs.dat @@ -0,0 +1,12 @@ +BEGIN options + DIGITS 7 + PRINT_INPUT +END options + +BEGIN continuous FILEOUT gwe_02.uze.obs.csv + mwe-1-temp TEMPERATURE 1 + mwe-1-stor STORAGE 1 + mwe-1-gwe1 UZE 1 + mwe-1-gwe2 UZE 2 + mwe-2-gwe1 UZE 3 +END continuous diff --git a/doc/mf6io/mf6ivar/examples/gwe-uze-example.dat b/doc/mf6io/mf6ivar/examples/gwe-uze-example.dat new file mode 100644 index 00000000000..e6e29baac43 --- /dev/null +++ b/doc/mf6io/mf6ivar/examples/gwe-uze-example.dat @@ -0,0 +1,24 @@ +BEGIN OPTIONS + AUXILIARY aux1 aux2 + BOUNDNAMES + PRINT_INPUT + PRINT_TEMPERATURE + PRINT_FLOWS + SAVE_FLOWS + TEMPERATURE FILEOUT gwe_02.uze.bin + BUDGET FILEOUT gwe_02.uze.bud + OBS6 FILEIN gwe_02.uze.obs +END OPTIONS + +BEGIN PACKAGEDATA +# L STRT aux1 aux2 bname + 1 0.0 99.0 999.0 MYUZFCELL1 + 2 0.0 99.0 999.0 MYUZFCELL2 + 3 0.0 99.0 999.0 MYUZFCELL3 +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 e30c47cf70b..05aa417e9ef 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -1448,6 +1448,40 @@ | GWE | SFE | PERIOD | AUXILIARY | KEYWORD | keyword for specifying auxiliary variable. | | GWE | SFE | 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 | SFE | 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 | UZE | 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 | UZE | 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 | UZE | 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 concentrations 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 | UZE | OPTIONS | BOUNDNAMES | KEYWORD | keyword to indicate that boundary names may be provided with the list of unsaturated zone flow cells. | +| GWE | UZE | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of unsaturated zone flow information will be written to the listing file immediately after it is read. | +| GWE | UZE | OPTIONS | PRINT_TEMPERATURE | KEYWORD | keyword to indicate that the list of UZF cell temperatures 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 temperatures are printed for the last time step of each stress period. | +| GWE | UZE | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of unsaturated zone 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 | UZE | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that unsaturated zone flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | +| GWE | UZE | OPTIONS | TEMPERATURE | KEYWORD | keyword to specify that record corresponds to temperature. | +| GWE | UZE | OPTIONS | TEMPFILE | STRING | name of the binary output file to write temperature information. | +| GWE | UZE | OPTIONS | BUDGET | KEYWORD | keyword to specify that record corresponds to the budget. | +| GWE | UZE | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. | +| GWE | UZE | OPTIONS | BUDGETFILE | STRING | name of the binary output file to write budget information. | +| GWE | UZE | OPTIONS | BUDGETCSV | KEYWORD | keyword to specify that record corresponds to the budget CSV. | +| GWE | UZE | 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 | UZE | OPTIONS | TS6 | KEYWORD | keyword to specify that record corresponds to a time-series file. | +| GWE | UZE | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | +| GWE | UZE | 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 | UZE | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. | +| GWE | UZE | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the UZE 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 UZE package. | +| GWE | UZE | PACKAGEDATA | UZFNO | INTEGER | integer value that defines the UZF cell number associated with the specified PACKAGEDATA data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. Unsaturated zone flow information must be specified for every UZF cell or the program will terminate with an error. The program also will terminate with an error if information for a UZF cell is specified more than once. | +| GWE | UZE | PACKAGEDATA | STRT | DOUBLE PRECISION | real value that defines the starting temperature for the unsaturated zone flow cell. | +| GWE | UZE | PACKAGEDATA | AUX | DOUBLE PRECISION (NAUX) | represents the values of the auxiliary variables for each unsaturated zone flow. The values of auxiliary variables must be present for each unsaturated zone flow. 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 | UZE | PACKAGEDATA | BOUNDNAME | STRING | name of the unsaturated zone flow 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 | UZE | 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 | UZE | PERIOD | UZFNO | INTEGER | integer value that defines the UZF cell number associated with the specified PERIOD data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. | +| GWE | UZE | PERIOD | UZESETTING | KEYSTRING | line of information that is parsed into a keyword and values. Keyword values that can be used to start the UZESETTING string include: STATUS, TEMPERATURE, INFILTRATION, UZET, and AUXILIARY. These settings are used to assign the temperature associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. | +| GWE | UZE | PERIOD | STATUS | STRING | keyword option to define UZF cell status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the UZF cell. If a UZF cell is inactive, then there will be no energy fluxes into or out of the UZF cell and the inactive value will be written for the UZF cell temperature. If a UZF cell is constant, then the temperature for the UZF cell will be fixed at the user specified value. | +| GWE | UZE | PERIOD | TEMPERATURE | STRING | real or character value that defines the temperature for the unsaturated zone flow cell. The specified TEMPERATURE is only applied if the unsaturated zone flow cell is a constant temperature cell. 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 | UZE | PERIOD | INFILTRATION | STRING | real or character value that defines the temperature of the infiltration $(^\circ C)$ for the UZF cell. 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 | UZE | PERIOD | UZET | STRING | real or character value that states what fraction of the simulated unsaturated zone evapotranspiration is associated with evaporation. The evaporative losses from the unsaturated zone moisture content will have an evaporative cooling effect on the GWE cell from which the evaporation originated. If this value is larger than 1, then it will be reset to 1. 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 | UZE | PERIOD | AUXILIARY | KEYWORD | keyword for specifying auxiliary variable. | +| GWE | UZE | 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 | UZE | 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 | FMI | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that FMI flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | | GWE | FMI | OPTIONS | FLOW_IMBALANCE_CORRECTION | KEYWORD | correct for an imbalance in flows by assuming that any residual flow error comes in or leaves at the temperature of the cell. When this option is activated, the GWE Model budget written to the listing file will contain two additional entries: FLOW-ERROR and FLOW-CORRECTION. These two entries will be equal but opposite in sign. The FLOW-CORRECTION term is a mass flow that is added to offset the error caused by an imprecise flow balance. If these terms are not relatively small, the flow model should be rerun with stricter convergence tolerances. | | GWE | FMI | PACKAGEDATA | FLOWTYPE | STRING | is the word GWFBUDGET, GWFHEAD, GWFMOVER or the name of an advanced GWF stress package. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. If an advanced GWF stress package name appears then the corresponding file must be the budget file saved by a LAK, SFR, MAW or UZF Package. | diff --git a/doc/mf6io/mf6ivar/mf6ivar.py b/doc/mf6io/mf6ivar/mf6ivar.py index 1170f4a7428..a968c2c0ff0 100644 --- a/doc/mf6io/mf6ivar/mf6ivar.py +++ b/doc/mf6io/mf6ivar/mf6ivar.py @@ -754,6 +754,7 @@ def write_appendix(texdir, allblocks): "gwe-oc", "gwe-ssm", "gwe-sfe", + "gwe-uze", "gwe-fmi", "utl-spc", "utl-spca", diff --git a/doc/mf6io/mf6ivar/tex/appendixA.tex b/doc/mf6io/mf6ivar/tex/appendixA.tex index cf214f3156f..b9581283c5e 100644 --- a/doc/mf6io/mf6ivar/tex/appendixA.tex +++ b/doc/mf6io/mf6ivar/tex/appendixA.tex @@ -302,6 +302,10 @@ GWE & SFE & PACKAGEDATA & yes \\ GWE & SFE & PERIOD & yes \\ \hline +GWE & UZE & OPTIONS & yes \\ +GWE & UZE & PACKAGEDATA & yes \\ +GWE & UZE & PERIOD & yes \\ +\hline GWE & FMI & OPTIONS & yes \\ GWE & FMI & PACKAGEDATA & yes \\ \hline diff --git a/doc/mf6io/mf6ivar/tex/gwe-uze-desc.tex b/doc/mf6io/mf6ivar/tex/gwe-uze-desc.tex new file mode 100644 index 00000000000..791b8dd9404 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-uze-desc.tex @@ -0,0 +1,91 @@ +% 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 concentrations 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 unsaturated zone flow cells. + +\item \texttt{PRINT\_INPUT}---keyword to indicate that the list of unsaturated zone flow information will be written to the listing file immediately after it is read. + +\item \texttt{PRINT\_TEMPERATURE}---keyword to indicate that the list of UZF cell temperatures 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 temperatures are printed for the last time step of each stress period. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that the list of unsaturated zone 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 unsaturated zone 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 UZE 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 UZE package. + +\end{description} +\item \textbf{Block: PACKAGEDATA} + +\begin{description} +\item \texttt{uzfno}---integer value that defines the UZF cell number associated with the specified PACKAGEDATA data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. Unsaturated zone flow information must be specified for every UZF cell or the program will terminate with an error. The program also will terminate with an error if information for a UZF cell is specified more than once. + +\item \texttt{strt}---real value that defines the starting temperature for the unsaturated zone flow cell. + +\item \textcolor{blue}{\texttt{aux}---represents the values of the auxiliary variables for each unsaturated zone flow. The values of auxiliary variables must be present for each unsaturated zone flow. 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 unsaturated zone flow 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{uzfno}---integer value that defines the UZF cell number associated with the specified PERIOD data on the line. UZFNO must be greater than zero and less than or equal to NUZFCELLS. + +\item \texttt{uzesetting}---line of information that is parsed into a keyword and values. Keyword values that can be used to start the UZESETTING string include: STATUS, TEMPERATURE, INFILTRATION, UZET, and AUXILIARY. These settings are used to assign the temperature associated with the corresponding flow terms. Temperatures cannot be specified for all flow terms. + +\begin{lstlisting}[style=blockdefinition] +STATUS +TEMPERATURE <@temperature@> +INFILTRATION <@infiltration@> +UZET <@uzet@> +AUXILIARY <@auxval@> +\end{lstlisting} + +\item \texttt{status}---keyword option to define UZF cell status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that temperature will be calculated for the UZF cell. If a UZF cell is inactive, then there will be no energy fluxes into or out of the UZF cell and the inactive value will be written for the UZF cell temperature. If a UZF cell is constant, then the temperature for the UZF cell will be fixed at the user specified value. + +\item \textcolor{blue}{\texttt{temperature}---real or character value that defines the temperature for the unsaturated zone flow cell. The specified TEMPERATURE is only applied if the unsaturated zone flow cell is a constant temperature cell. 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{infiltration}---real or character value that defines the temperature of the infiltration $(^\circ C)$ for the UZF cell. 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{uzet}---real or character value that states what fraction of the simulated unsaturated zone evapotranspiration is associated with evaporation. The evaporative losses from the unsaturated zone moisture content will have an evaporative cooling effect on the GWE cell from which the evaporation originated. If this value is larger than 1, then it will be reset to 1. 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-uze-options.dat b/doc/mf6io/mf6ivar/tex/gwe-uze-options.dat new file mode 100644 index 00000000000..ef83c7fa718 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-uze-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-uze-packagedata.dat b/doc/mf6io/mf6ivar/tex/gwe-uze-packagedata.dat new file mode 100644 index 00000000000..b6b04c298fe --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-uze-packagedata.dat @@ -0,0 +1,5 @@ +BEGIN PACKAGEDATA + [<@aux(naux)@>] [] + [<@aux(naux)@>] [] + ... +END PACKAGEDATA diff --git a/doc/mf6io/mf6ivar/tex/gwe-uze-period.dat b/doc/mf6io/mf6ivar/tex/gwe-uze-period.dat new file mode 100644 index 00000000000..b06edb21028 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/gwe-uze-period.dat @@ -0,0 +1,5 @@ +BEGIN PERIOD + + + ... +END PERIOD diff --git a/make/makefile b/make/makefile index 07f580b3900..6b71b70c792 100644 --- a/make/makefile +++ b/make/makefile @@ -279,6 +279,7 @@ $(OBJDIR)/gwf3buy8.o \ $(OBJDIR)/GhostNode.o \ $(OBJDIR)/gwf3evt8.o \ $(OBJDIR)/gwf3chd8.o \ +$(OBJDIR)/gwe1uze1.o \ $(OBJDIR)/gwe1sfe1.o \ $(OBJDIR)/gwe1mwe1.o \ $(OBJDIR)/gwe1lke1.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index ad8ea63a9d0..32f92daee02 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -154,7 +154,8 @@ - + + diff --git a/src/Model/GroundWaterEnergy/gwe1.f90 b/src/Model/GroundWaterEnergy/gwe1.f90 index dd5536dba78..fafde6d3101 100644 --- a/src/Model/GroundWaterEnergy/gwe1.f90 +++ b/src/Model/GroundWaterEnergy/gwe1.f90 @@ -765,7 +765,7 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & use GweLkeModule, only: lke_create use GweSfeModule, only: sfe_create use GweMweModule, only: mwe_create - !use GweUzeModule, only: uze_create + use GweUzeModule, only: uze_create use ApiModule, only: api_create ! -- dummy class(GweModelType) :: this @@ -802,10 +802,10 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, & 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, & - ! this%gwecommon) + case ('UZE6') + call uze_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & + pakname, this%fmi, this%eqnsclfac, this%gwecommon, & + this%depvartype, this%depvarunit, this%depvarunitabbrev) !case('API6') ! call api_create(packobj, ipakid, ipaknum, inunit, iout, this%name, & ! pakname) diff --git a/src/Model/GroundWaterEnergy/gwe1uze1.f90 b/src/Model/GroundWaterEnergy/gwe1uze1.f90 new file mode 100644 index 00000000000..64231abc2ff --- /dev/null +++ b/src/Model/GroundWaterEnergy/gwe1uze1.f90 @@ -0,0 +1,1395 @@ +! -- Unsaturated Zone Flow Energy Transport Module +! -- todo: save the uze temperature into the uze aux variable? +! -- todo: calculate the uzf DENSE aux variable using temperature? +! -- todo: GWD and GWD-TO-MVR do not seem to be included; prob in UZF? +! +! UZF flows (flowbudptr) index var UZE term Transport Type +!--------------------------------------------------------------------------------- + +! -- terms from UZF that will be handled by parent APT Package +! FLOW-JA-FACE idxbudfjf FLOW-JA-FACE cv2cv +! GWF (aux FLOW-AREA) idxbudgwf GWF uzf2gwf +! STORAGE (aux VOLUME) idxbudsto none used for water volumes +! FROM-MVR idxbudfmvr FROM-MVR q * cext = this%qfrommvr(:) +! AUXILIARY none none none +! none none STORAGE (aux MASS) +! none none AUXILIARY none + +! -- terms from UZF that need to be handled here +! INFILTRATION idxbudinfl INFILTRATION q < 0: q * cwell, else q * cuser +! REJ-INF idxbudrinf REJ-INF q * cuze +! UZET idxbuduzet UZET q * cet +! REJ-INF-TO-MVR idxbudritm REJ-INF-TO-MVR q * cinfil? +! THERMAL-EQUIL idxbudtheq THERMAL-EQUIL residual + +! -- terms from UZF that should be skipped + +module GweUzeModule + + use KindModule, only: DP, I4B + use ConstantsModule, only: DZERO, DONE, LINELENGTH + use SimModule, only: store_error + use BndModule, only: BndType, GetBndFromList + use TspFmiModule, only: TspFmiType + use UzfModule, only: UzfType + use ObserveModule, only: ObserveType + use TspAptModule, only: TspAptType, apt_process_obsID, & + apt_process_obsID12 + use GweInputDataModule, only: GweInputDataType + use MatrixBaseModule + + implicit none + + public uze_create + + character(len=*), parameter :: ftype = 'UZE' + character(len=*), parameter :: flowtype = 'UZF' + character(len=16) :: text = ' UZE' + + type, extends(TspAptType) :: GweUzeType + + type(GweInputDataType), pointer :: gwecommon => null() !< pointer to shared gwe data used by multiple packages but set in mst + + integer(I4B), pointer :: idxbudinfl => null() !< index of uzf infiltration terms in flowbudptr + integer(I4B), pointer :: idxbudrinf => null() !< index of rejected infiltration terms in flowbudptr + integer(I4B), pointer :: idxbuduzet => null() !< index of unsat et terms in flowbudptr + integer(I4B), pointer :: idxbudritm => null() !< index of rej infil to mover rate to mover terms in flowbudptr + integer(I4B), pointer :: idxbudtheq => null() !< index of thermal equilibration terms in flowbudptr + + real(DP), dimension(:), pointer, contiguous :: tempinfl => null() !< infiltration temperature + real(DP), dimension(:), pointer, contiguous :: tempuzet => null() !< unsat et temperature + + contains + + procedure :: bnd_da => uze_da + procedure :: allocate_scalars + procedure :: apt_allocate_arrays => uze_allocate_arrays + procedure :: find_apt_package => find_uze_package + procedure :: apt_fc_expanded => uze_fc_expanded + procedure :: pak_solve => uze_solve + procedure :: pak_get_nbudterms => uze_get_nbudterms + procedure :: pak_setup_budobj => uze_setup_budobj + procedure :: pak_fill_budobj => uze_fill_budobj + procedure :: uze_infl_term + procedure :: uze_rinf_term + procedure :: uze_uzet_term + procedure :: uze_ritm_term + procedure :: uze_theq_term + procedure :: pak_df_obs => uze_df_obs + procedure :: pak_rp_obs => uze_rp_obs + procedure :: pak_bd_obs => uze_bd_obs + procedure :: pak_set_stressperiod => uze_set_stressperiod + procedure :: bnd_ac => uze_ac + procedure :: bnd_mc => uze_mc + + end type GweUzeType + +contains + + !> @breif Create a new UZE package + !< + subroutine uze_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(GweUzeType), pointer :: uzeobj + ! + ! -- Allocate the object and assign values to object variables + allocate (uzeobj) + packobj => uzeobj + ! + ! -- Create name and memory path + call packobj%set_names(ibcnum, namemodel, pakname, ftype) + packobj%text = text + ! + ! -- Allocate scalars + call uzeobj%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 GwfGwt exchange is + ! created, it sets fmi%bndlist so that the GWT model has access to all + ! the flow packages + uzeobj%fmi => fmi + ! + ! -- Store pointer to governing equation scale factor + uzeobj%eqnsclfac => eqnsclfac + ! + ! -- Store pointer to shared data module for accessing cpw, rhow + ! for the budget calculations, and for accessing the latent heat of + ! vaporization + uzeobj%gwecommon => gwecommon + ! + ! -- Set labels that will be used in generalized APT class + uzeobj%depvartype = dvt + uzeobj%depvarunit = dvu + uzeobj%depvarunitabbrev = dvua + ! + ! -- Return + return + end subroutine uze_create + + !> @brief Find corresponding uze package + !< + subroutine find_uze_package(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweUzeType) :: 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 (UzfType) + 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 ('INFILTRATION') + this%idxbudinfl = ip + this%idxbudssm(ip) = 0 + case ('REJ-INF') + this%idxbudrinf = ip + this%idxbudssm(ip) = 0 + case ('UZET') + this%idxbuduzet = ip + this%idxbudssm(ip) = 0 + case ('REJ-INF-TO-MVR') + this%idxbudritm = 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 temperatures + ! 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' + ! + ! -- Thermal equilibration term + this%idxbudtheq = this%flowbudptr%nbudterm + 1 + ! + ! -- Return + return + end subroutine find_uze_package + + !> @brief Add package connection to matrix. + !! + !! Overrides apt_ac to fold the UZE heat balance terms into the row + !! corresponding to the host cell and enforce thermal equilibrium between + !! UZE and the GWE cell. + !< + subroutine uze_ac(this, moffset, sparse) + use MemoryManagerModule, only: mem_setptr + use SparseModule, only: sparsematrix + ! -- dummy + class(GweUzeType), intent(inout) :: this + integer(I4B), intent(in) :: moffset !< current models starting position in global matrix + type(sparsematrix), intent(inout) :: sparse + ! -- local + integer(I4B) :: i, ii + integer(I4B) :: n !< index of a uze object within the complete list of uze objects + integer(I4B) :: jj !< + integer(I4B) :: nglo !< index of uze object in global matrix + integer(I4B) :: jglo !< host cell's position in global matrix for a uze object + integer(I4B) :: idxn !< used for cross-check + integer(I4B) :: idxjj !< used for cross-check + integer(I4B) :: idxnglo !< used for cross-check + integer(I4B) :: idxjglo !< used for cross-check + ! + ! -- Add package rows to sparse + if (this%imatrows /= 0) then + ! + ! -- Diagonal on the row assoc. with the uze feature + do n = 1, this%ncv + nglo = moffset + this%dis%nodes + this%ioffset + n + call sparse%addconnection(nglo, nglo, 1) + end do + ! + ! -- Add uze-to-gwe connections. For uze, this particular do loop + ! is the same as its counterpart in apt_ac. + ! nlist: number of gwe cells with a connection to at least one uze object + do i = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + n = this%flowbudptr%budterm(this%idxbudgwf)%id1(i) !< uze object position within uze object list + jj = this%flowbudptr%budterm(this%idxbudgwf)%id2(i) !< position of gwe cell to which uze feature is connected + nglo = moffset + this%dis%nodes + this%ioffset + n !< uze feature position + jglo = moffset + jj !< gwe cell position + call sparse%addconnection(nglo, jglo, 1) + call sparse%addconnection(jglo, nglo, 1) + end do + ! + ! -- For uze, add feature-to-feature connections (i.e., + ! vertically stacked UZ objects) to row corresponding + ! to the host cell. Terms added to the row associated with host + ! cell are added to columns associated with other uze features. + ! This approach deviates from the approach taken in apt_ac. + if (this%idxbudfjf /= 0) then + do i = 1, this%flowbudptr%budterm(this%idxbudfjf)%maxlist + n = this%flowbudptr%budterm(this%idxbudfjf)%id1(i) !< position of currently considered uze feature + jj = this%flowbudptr%budterm(this%idxbudfjf)%id2(i) !< position of connected uze feature + nglo = moffset + this%dis%nodes + this%ioffset + n !< global position of currently considered uze feature + jglo = moffset + this%dis%nodes + this%ioffset + jj !< global position of connected uze feature + ! -- If connected uze feature is upstream, find cell that hosts currently + ! considered uze feature and add connection to that cell's row + do ii = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist !< uze object id among uze objects + idxn = this%flowbudptr%budterm(this%idxbudgwf)%id1(ii) !< uze object position within uze object list + idxjj = this%flowbudptr%budterm(this%idxbudgwf)%id2(ii) !< position of gwe cell to which uze feature is connected + idxnglo = moffset + this%dis%nodes + this%ioffset + idxn !< uze feature global position + idxjglo = moffset + idxjj !< gwe cell global position + if (nglo == idxnglo) exit + end do + call sparse%addconnection(idxjglo, jglo, 1) + end do + end if + end if + ! + ! -- Return + return + end subroutine uze_ac + + !> @brief Map package connection to matrix + !< + subroutine uze_mc(this, moffset, matrix_sln) + use SparseModule, only: sparsematrix + ! -- dummy + class(GweUzeType), intent(inout) :: this + integer(I4B), intent(in) :: moffset + class(MatrixBaseType), pointer :: matrix_sln + ! -- local + integer(I4B) :: n, j, iglo, jglo + integer(I4B) :: idxn, idxj, idxiglo, idxjglo + integer(I4B) :: ipos, idxpos + ! + ! -- Allocate memory for index arrays + call this%apt_allocate_index_arrays() + ! + ! -- Store index positions + if (this%imatrows /= 0) then + ! + ! -- Find the position of each connection in the global ia, ja structure + ! and store them in idxglo. idxglo allows this model to insert or + ! retrieve values into or from the global A matrix + ! apt rows + ! + ! -- Feature diagonal in global matrix + do n = 1, this%ncv + iglo = moffset + this%dis%nodes + this%ioffset + n + this%idxpakdiag(n) = matrix_sln%get_position_diag(iglo) + end do + ! + ! -- Cell to feature connection in global matrix + do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos) !< feature number + j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos) !< cell number + iglo = moffset + this%dis%nodes + this%ioffset + n !< feature row index + jglo = j + moffset !< cell row index + ! -- Note that this is where idxlocnode is set for uze; it is set + ! to the host cell local row index rather than the feature local + ! row index + this%idxlocnode(n) = j ! kluge note: do we want to introduce a new array instead of co-opting idxlocnode??? + ! -- For connection ipos in list of feature-cell connections, + ! global positions of feature-row diagonal and off-diagonal + ! corresponding to the cell + this%idxdglo(ipos) = matrix_sln%get_position_diag(iglo) + this%idxoffdglo(ipos) = matrix_sln%get_position(iglo, jglo) + end do + ! + ! -- Feature to cell connection in global matrix + do ipos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + n = this%flowbudptr%budterm(this%idxbudgwf)%id1(ipos) !< feature number + j = this%flowbudptr%budterm(this%idxbudgwf)%id2(ipos) !< cell number + iglo = j + moffset !< cell row index + jglo = moffset + this%dis%nodes + this%ioffset + n !< feature row index + ! -- For connection ipos in list of feature-cell connections, + ! global positions of cell-row diagonal and off-diagonal + ! corresponding to the feature + this%idxsymdglo(ipos) = matrix_sln%get_position_diag(iglo) + this%idxsymoffdglo(ipos) = matrix_sln%get_position(iglo, jglo) + end do + ! + ! -- Feature to feature connection in global matrix + if (this%idxbudfjf /= 0) then + do ipos = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist + n = this%flowbudptr%budterm(this%idxbudfjf)%id1(ipos) !< number of currently considered uze feature + j = this%flowbudptr%budterm(this%idxbudfjf)%id2(ipos) !< number of connected uze feature + iglo = moffset + this%dis%nodes + this%ioffset + n !< global position of currently considered uze feature + jglo = moffset + this%dis%nodes + this%ioffset + j !< global position of connected uze feature + ! -- If connected uze feature is upstream, find cell that hosts currently + ! considered uze feature and map connection to that cell's row + do idxpos = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + idxn = this%flowbudptr%budterm(this%idxbudgwf)%id1(idxpos) !< feature number + idxj = this%flowbudptr%budterm(this%idxbudgwf)%id2(idxpos) !< cell number) + idxjglo = moffset + this%dis%nodes + this%ioffset + idxn !< feature row index + idxiglo = moffset + idxj !< cell row index + if (idxjglo == iglo) exit + end do + ! -- For connection ipos in list of feature-feature connections, + ! global positions of host-cell-row entries corresponding to + ! (in the same columns as) the feature-id1-row diagonal and the + ! feature-id1-row off-diagonal corresponding to feature id2 + this%idxfjfdglo(ipos) = matrix_sln%get_position_diag(idxiglo) + this%idxfjfoffdglo(ipos) = matrix_sln%get_position(idxiglo, jglo) + end do + end if + end if + ! + ! -- Return + return + end subroutine uze_mc + + !> @brief Add matrix terms related to UZE + !! + !! This will be called from TspAptType%apt_fc_expanded() + !! in order to add matrix terms specifically for this package + !< + subroutine uze_fc_expanded(this, rhs, ia, idxglo, matrix_sln) + ! -- dummy + class(GweUzeType) :: 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) :: ipossymoffd + real(DP) :: cold + real(DP) :: qbnd + real(DP) :: omega + real(DP) :: rrate + real(DP) :: rhsval + real(DP) :: hcofval + ! + ! -- Add infiltration contribution + ! uze does not put feature balance coefficients in the row + ! associated with the feature. Instead, these coefficients are + ! moved into the row associated with cell hosting the uze feature + if (this%idxbudinfl /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudinfl)%nlist + call this%uze_infl_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) !< for uze idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add rejected infiltration contribution + if (this%idxbudrinf /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrinf)%nlist + call this%uze_rinf_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add unsaturated et contribution + if (this%idxbuduzet /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbuduzet)%nlist + call this%uze_uzet_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add rejected infiltration to mover contribution + if (this%idxbudritm /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudritm)%nlist + call this%uze_ritm_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) ! for uze idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- For UZE, content of apt_fc_expanded placed here as the approach is to + ! completely override apt_fc_expanded() with what follows + ! + ! -- Mass (or energy) storage in features + do n = 1, this%ncv + cold = this%xoldpak(n) + iloc = this%idxlocnode(n) ! for uze idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(n) + call this%apt_stor_term(n, n1, n2, rrate, rhsval, hcofval) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + ! + ! -- Add to mover contribution + if (this%idxbudtmvr /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudtmvr)%nlist + call this%apt_tmvr_term(j, n1, n2, rrate, rhsval, hcofval) + iloc = this%idxlocnode(n1) ! for uze, idxlocnode stores the host cell local row index + ipossymoffd = this%idxsymoffdglo(j) + call matrix_sln%add_value_pos(ipossymoffd, hcofval) + rhs(iloc) = rhs(iloc) + rhsval + end do + end if + ! + ! -- Add from mover contribution + if (this%idxbudfmvr /= 0) then + do n = 1, this%ncv + rhsval = this%qmfrommvr(n) ! kluge note: presumably already in terms of energy + iloc = this%idxlocnode(n) ! for uze idxlocnode stores the host cell local row index + rhs(iloc) = rhs(iloc) - rhsval + end do + end if + ! + ! -- Go through each apt-gwf connection + do j = 1, this%flowbudptr%budterm(this%idxbudgwf)%nlist + ! + ! -- Set n to feature number and process if active feature + n = this%flowbudptr%budterm(this%idxbudgwf)%id1(j) + if (this%iboundpak(n) /= 0) then + ! + ! -- This code altered from its counterpart appearing in apt; this equates + ! uze temperature to cell temperature using the feature's row + iposd = this%idxdglo(j) + iposoffd = this%idxoffdglo(j) + call matrix_sln%add_value_pos(iposd, DONE) + call matrix_sln%add_value_pos(iposoffd, -DONE) + end if + end do + ! + ! -- Go through each apt-apt connection + if (this%idxbudfjf /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudfjf)%nlist + n1 = this%flowbudptr%budterm(this%idxbudfjf)%id1(j) + n2 = this%flowbudptr%budterm(this%idxbudfjf)%id2(j) + qbnd = this%flowbudptr%budterm(this%idxbudfjf)%flow(j) + if (qbnd <= DZERO) then + omega = DONE + else + omega = DZERO + end if + iposd = this%idxfjfdglo(j) !< position of feature-id1 column in feature id1's host-cell row + iposoffd = this%idxfjfoffdglo(j) !< position of feature-id2 column in feature id1's host-cell row + call matrix_sln%add_value_pos(iposd, omega * qbnd * this%eqnsclfac) + call matrix_sln%add_value_pos(iposoffd, & + (DONE - omega) * qbnd * this%eqnsclfac) + end do + end if + ! + ! -- Return + return + end subroutine uze_fc_expanded + + !> @brief Explicit solve + !! + !! There should be no explicit solve for uze. However, if there were, then + !! this subroutine would add terms specific to the unsaturated zone to the + !! explicit unsaturated-zone solve + subroutine uze_solve(this) + ! -- dummy + class(GweUzeType) :: this + ! -- local + integer(I4B) :: j + integer(I4B) :: n1, n2 + real(DP) :: rrate + ! + ! -- Add infiltration contribution + if (this%idxbudinfl /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudinfl)%nlist + call this%uze_infl_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add rejected infiltration contribution + if (this%idxbudrinf /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudrinf)%nlist + call this%uze_rinf_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add unsaturated et contribution + if (this%idxbuduzet /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbuduzet)%nlist + call this%uze_uzet_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Add rejected infiltration to mover contribution + if (this%idxbudritm /= 0) then + do j = 1, this%flowbudptr%budterm(this%idxbudritm)%nlist + call this%uze_ritm_term(j, n1, n2, rrate) + this%dbuff(n1) = this%dbuff(n1) + rrate + end do + end if + ! + ! -- Return + return + end subroutine uze_solve + + !> @brief Return the number of UZE-specific budget terms + !! + !! Function to return the number of budget terms just for this package. + !! This overrides function in parent. + !< + function uze_get_nbudterms(this) result(nbudterms) + ! -- dummy + class(GweUzeType) :: this + ! -- Return + integer(I4B) :: nbudterms + ! + ! -- Number of budget terms is 5 + nbudterms = 0 + if (this%idxbudinfl /= 0) nbudterms = nbudterms + 1 + if (this%idxbudrinf /= 0) nbudterms = nbudterms + 1 + if (this%idxbuduzet /= 0) nbudterms = nbudterms + 1 + if (this%idxbudritm /= 0) nbudterms = nbudterms + 1 + if (this%idxbudtheq /= 0) nbudterms = nbudterms + 1 + ! + ! -- Return + return + end function uze_get_nbudterms + + !> @brief Setup budget object + !! + !! Set up the budget object that stores all the unsaturated-zone + !! flows + !< + subroutine uze_setup_budobj(this, idx) + ! -- modules + use ConstantsModule, only: LENBUDTXT + ! -- dummy + class(GweUzeType) :: this + integer(I4B), intent(inout) :: idx + ! -- local + integer(I4B) :: maxlist, naux + character(len=LENBUDTXT) :: text + ! + ! -- Infiltration + text = ' INFILTRATION' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudinfl)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + ! + ! -- Rejected infiltration (Hortonian flow) + if (this%idxbudrinf /= 0) then + text = ' REJ-INF' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudrinf)%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 + ! + ! -- Evapotranspiration from the unsaturated zone + if (this%idxbuduzet /= 0) then + text = ' UZET' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbuduzet)%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 + ! + ! -- Rejected infiltration that is subsequently transferred by MVR + if (this%idxbudritm /= 0) then + text = ' INF-REJ-TO-MVR' + idx = idx + 1 + maxlist = this%flowbudptr%budterm(this%idxbudritm)%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 + ! + ! -- Energy transferred to solid phase by the thermal equilibrium assumption + text = ' THERMAL-EQUIL' + idx = idx + 1 + ! -- use dimension of GWF term + maxlist = this%flowbudptr%budterm(this%idxbudgwf)%maxlist + naux = 0 + call this%budobj%budterm(idx)%initialize(text, & + this%name_model, & + this%packName, & + this%name_model, & + this%packName, & + maxlist, .false., .false., & + naux) + ! + ! -- Return + return + end subroutine uze_setup_budobj + + !> @brief Fill UZE budget object + !! + !! Copy flow terms into this%budobj + !< + subroutine uze_fill_budobj(this, idx, x, flowja, ccratin, ccratout) + ! -- modules + ! -- dummy + class(GweUzeType) :: 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, indx + integer(I4B) :: nlist, nlen + integer(I4B) :: igwfnode + integer(I4B) :: idiag + real(DP) :: q + real(DP), dimension(:), allocatable :: budresid + ! + allocate (budresid(this%ncv)) + do n1 = 1, this%ncv + budresid(n1) = DZERO + end do + ! + indx = 0 + ! + ! -- FLOW JA FACE into budresid + nlen = 0 + if (this%idxbudfjf /= 0) then + nlen = this%flowbudptr%budterm(this%idxbudfjf)%maxlist + end if + if (nlen > 0) then + indx = indx + 1 + nlist = this%budobj%budterm(indx)%nlist + do j = 1, nlist + n1 = this%budobj%budterm(indx)%id1(j) + n2 = this%budobj%budterm(indx)%id2(j) + if (n1 < n2) then + q = this%budobj%budterm(indx)%flow(j) + budresid(n1) = budresid(n1) + q + budresid(n2) = budresid(n2) - q + end if + end do + end if + ! + ! -- GWF (LEAKAGE) into budresid + indx = indx + 1 + nlist = this%budobj%budterm(indx)%nlist + do j = 1, nlist + n1 = this%budobj%budterm(indx)%id1(j) + q = this%budobj%budterm(indx)%flow(j) + budresid(n1) = budresid(n1) + q + end do + ! + ! -- Skip individual package terms + indx = this%idxlastpak + ! + ! -- STORAGE into budresid + indx = indx + 1 + do n1 = 1, this%ncv + q = this%budobj%budterm(indx)%flow(n1) + budresid(n1) = budresid(n1) + q + end do + ! + ! -- TO MOVER into budresid + if (this%idxbudtmvr /= 0) then + indx = indx + 1 + nlist = this%budobj%budterm(indx)%nlist + do j = 1, nlist + n1 = this%budobj%budterm(indx)%id1(j) + q = this%budobj%budterm(indx)%flow(j) + budresid(n1) = budresid(n1) + q + end do + end if + ! + ! -- FROM MOVER into budresid + if (this%idxbudfmvr /= 0) then + indx = indx + 1 + nlist = this%budobj%budterm(indx)%nlist + do j = 1, nlist + n1 = this%budobj%budterm(indx)%id1(j) + q = this%budobj%budterm(indx)%flow(j) + budresid(n1) = budresid(n1) + q + end do + end if + ! + ! -- CONSTANT FLOW into budresid + indx = indx + 1 + do n1 = 1, this%ncv + q = this%budobj%budterm(indx)%flow(n1) + budresid(n1) = budresid(n1) + q + end do + ! + ! -- AUXILIARY VARIABLES into budresid + ! -- (No flows associated with these) + ! + ! -- Individual package terms processed last + ! + ! -- Infiltration + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudinfl)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%uze_infl_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + budresid(n1) = budresid(n1) + q + end do + ! + ! -- Rej-Inf + if (this%idxbudrinf /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudrinf)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%uze_rinf_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + budresid(n1) = budresid(n1) + q + end do + end if + ! + ! -- UZET + if (this%idxbuduzet /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbuduzet)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%uze_uzet_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + budresid(n1) = budresid(n1) + q + end do + end if + ! + ! -- Rej-Inf-To-MVR + if (this%idxbudritm /= 0) then + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudritm)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + call this%uze_ritm_term(j, n1, n2, q) + call this%budobj%budterm(idx)%update_term(n1, n2, q) + call this%apt_accumulate_ccterm(n1, q, ccratin, ccratout) + budresid(n1) = budresid(n1) + q + end do + end if + ! + ! -- Thermal-Equil + ! -- processed last because it is calculated from the residual + idx = idx + 1 + nlist = this%flowbudptr%budterm(this%idxbudgwf)%nlist + call this%budobj%budterm(idx)%reset(nlist) + do j = 1, nlist + n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(j) + igwfnode = this%flowbudptr%budterm(this%idxbudgwf)%id2(j) + q = -budresid(n1) + call this%uze_theq_term(j, n1, igwfnode, q) + 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(n1) = this%simvals(n1) - q + idiag = this%dis%con%ia(igwfnode) + flowja(idiag) = flowja(idiag) - q + end if + end do + ! + deallocate (budresid) + ! + ! -- Return + return + end subroutine uze_fill_budobj + + !> @brief Allocate scalars + !! + !! Allocate scalars specific to UZE package + !< + subroutine allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweUzeType) :: this + ! -- local + ! + ! -- Allocate scalars in TspAptType + call this%TspAptType%allocate_scalars() + ! + ! -- Allocate + call mem_allocate(this%idxbudinfl, 'IDXBUDINFL', this%memoryPath) + call mem_allocate(this%idxbudrinf, 'IDXBUDRINF', this%memoryPath) + call mem_allocate(this%idxbuduzet, 'IDXBUDUZET', this%memoryPath) + call mem_allocate(this%idxbudritm, 'IDXBUDRITM', this%memoryPath) + call mem_allocate(this%idxbudtheq, 'IDXBUDTHEQ', this%memoryPath) + ! + ! -- Initialize + this%idxbudinfl = 0 + this%idxbudrinf = 0 + this%idxbuduzet = 0 + this%idxbudritm = 0 + this%idxbudtheq = 0 + ! + ! -- Return + return + end subroutine allocate_scalars + + !> @brief Allocate arrays + !! + !! Allocate arrays used by the UZE package + !< + subroutine uze_allocate_arrays(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(GweUzeType), intent(inout) :: this + ! -- local + integer(I4B) :: n + ! + ! -- Time series + call mem_allocate(this%tempinfl, this%ncv, 'TEMPINFL', this%memoryPath) + call mem_allocate(this%tempuzet, this%ncv, 'TEMPUZET', this%memoryPath) + ! + ! -- Call standard TspAptType allocate arrays + call this%TspAptType%apt_allocate_arrays() + ! + ! -- Initialize + do n = 1, this%ncv + this%tempinfl(n) = DZERO + this%tempuzet(n) = DZERO + end do + ! + ! -- Return + return + end subroutine uze_allocate_arrays + + !> @brief Deallocate memory + !< + subroutine uze_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(GweUzeType) :: this + ! + ! -- Deallocate scalars + call mem_deallocate(this%idxbudinfl) + call mem_deallocate(this%idxbudrinf) + call mem_deallocate(this%idxbuduzet) + call mem_deallocate(this%idxbudritm) + call mem_deallocate(this%idxbudtheq) + ! + ! -- Deallocate time series + call mem_deallocate(this%tempinfl) + call mem_deallocate(this%tempuzet) + ! + ! -- Deallocate scalars in TspAptType + call this%TspAptType%bnd_da() + ! + ! -- Return + return + end subroutine uze_da + + !> @brief Infiltration term + !! + !! Accounts for energy added to the subsurface via infiltration, for example, + !! energy entering the model domain via rainfall or irrigation. + !< + subroutine uze_infl_term(this, ientry, n1, n2, rrate, & + rhsval, hcofval) + ! -- dummy + class(GweUzeType) :: 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%idxbudinfl)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudinfl)%id2(ientry) + ! + ! -- Note that qbnd is negative for negative infiltration + qbnd = this%flowbudptr%budterm(this%idxbudinfl)%flow(ientry) + if (qbnd < DZERO) then + ctmp = this%xnewpak(n1) + h = qbnd + r = DZERO + else + ctmp = this%tempinfl(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 uze_infl_term + + !> @brief Rejected infiltration term + !! + !! Accounts for energy that is added to the model from specifying an + !! infiltration rate and temperature, but is subsequently removed from + !! the model as that portion of the infiltration that is rejected (and + !! NOT transferred to another advanced package via the MVR/MVT packages). + !< + subroutine uze_rinf_term(this, ientry, n1, n2, rrate, & + rhsval, hcofval) + ! -- dummy + class(GweUzeType) :: 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%idxbudrinf)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudrinf)%id2(ientry) + qbnd = this%flowbudptr%budterm(this%idxbudrinf)%flow(ientry) + ctmp = this%tempinfl(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 uze_rinf_term + + !> @brief Evapotranspiration from the unsaturated-zone term + !! + !! Accounts for thermal cooling in the unsaturated zone as a result of + !! evapotranspiration from the unsaturated zone. Amount of water converted + !! to vapor phase (UZET) determined by GWF model + !< + subroutine uze_uzet_term(this, ientry, n1, n2, rrate, rhsval, hcofval) + ! -- dummy + class(GweUzeType) :: 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) :: omega + ! + n1 = this%flowbudptr%budterm(this%idxbuduzet)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbuduzet)%id2(ientry) + ! -- Note that qbnd is negative for uzet + qbnd = this%flowbudptr%budterm(this%idxbuduzet)%flow(ientry) + ctmp = this%tempuzet(n1) + if (this%xnewpak(n1) < ctmp) then + omega = DONE + else + omega = DZERO + end if + if (present(rrate)) & + rrate = (omega * qbnd * this%xnewpak(n1) + & + (DONE - omega) * qbnd * ctmp) * this%eqnsclfac + if (present(rhsval)) rhsval = -(DONE - omega) * qbnd * ctmp * this%eqnsclfac + if (present(hcofval)) hcofval = omega * qbnd * this%eqnsclfac + ! + ! -- Return + return + end subroutine uze_uzet_term + + !> @brief Rejected infiltration to MVR/MVT term + !! + !! Accounts for energy that is added to the model from specifying an + !! infiltration rate and temperature, but does not infiltrate into the + !! subsurface. This subroutine is called when the rejected infiltration + !! is transferred to another advanced package via the MVR/MVT packages. + !< + subroutine uze_ritm_term(this, ientry, n1, n2, rrate, & + rhsval, hcofval) + ! -- dummy + class(GweUzeType) :: 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%idxbudritm)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudritm)%id2(ientry) + qbnd = this%flowbudptr%budterm(this%idxbudritm)%flow(ientry) + ctmp = this%tempinfl(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 uze_ritm_term + + !> @brief Heat transferred through thermal equilibrium with the solid phase + !! + !! Accounts for the transfer of energy from the liquid phase to the solid + !! phase as a result of the instantaneous thermal equilibrium assumption. + !< + subroutine uze_theq_term(this, ientry, n1, n2, rrate) + ! -- modules + use ConstantsModule, only: LENBUDTXT + ! -- dummy + class(GweUzeType) :: this + integer(I4B), intent(in) :: ientry + integer(I4B), intent(inout) :: n1 + integer(I4B), intent(inout) :: n2 + real(DP), intent(inout) :: rrate + ! -- local + real(DP) :: r + integer(I4B) :: i + character(len=LENBUDTXT) :: flowtype + ! + r = DZERO + n1 = this%flowbudptr%budterm(this%idxbudgwf)%id1(ientry) + n2 = this%flowbudptr%budterm(this%idxbudgwf)%id2(ientry) + if (this%iboundpak(n1) /= 0) then + do i = 1, this%budobj%nbudterm + flowtype = this%budobj%budterm(i)%flowtype + select case (trim(adjustl(flowtype))) + case ('THERMAL-EQUIL') + ! -- Skip + continue + case default + r = r - this%budobj%budterm(i)%flow(ientry) + end select + end do + end if + rrate = r + ! + ! -- Return + return + end subroutine uze_theq_term + + !> @brief Define UZE Observation + !! + !! This subroutine: + !! - Stores observation types supported by the parent APT package. + !! - Overrides BndType%bnd_df_obs + !< + subroutine uze_df_obs(this) + ! -- dummy + class(GweUzeType) :: 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 + ! + ! -- Store obs type and assign procedure pointer + ! for flow between uze cells. + call this%obs%StoreObsType('flow-ja-face', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID12 + ! + ! -- 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 uze + !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: uze + call this%obs%StoreObsType('uze', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('infiltration', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('rej-inf', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('uzet', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Store obs type and assign procedure pointer + ! for observation type. + call this%obs%StoreObsType('rej-inf-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('thermal-equil', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => apt_process_obsID + ! + ! -- Return + return + end subroutine uze_df_obs + + !> @brief Process package specific obs + !! + !! Method to process specific observations for this package. + !< + subroutine uze_rp_obs(this, obsrv, found) + ! -- dummy + class(GweUzeType), 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 ('INFILTRATION') + call this%rp_obs_byfeature(obsrv) + case ('REJ-INF') + call this%rp_obs_byfeature(obsrv) + case ('UZET') + call this%rp_obs_byfeature(obsrv) + case ('REJ-INF-TO-MVR') + call this%rp_obs_byfeature(obsrv) + case ('THERMAL-EQUIL') + call this%rp_obs_byfeature(obsrv) + case default + found = .false. + end select + ! + return + end subroutine uze_rp_obs + + !> @brief Calculate observation value and pass it back to APT + !< + subroutine uze_bd_obs(this, obstypeid, jj, v, found) + ! -- dummy + class(GweUzeType), 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 ('INFILTRATION') + if (this%iboundpak(jj) /= 0 .and. this%idxbudinfl > 0) then + call this%uze_infl_term(jj, n1, n2, v) + end if + case ('REJ-INF') + if (this%iboundpak(jj) /= 0 .and. this%idxbudrinf > 0) then + call this%uze_rinf_term(jj, n1, n2, v) + end if + case ('UZET') + if (this%iboundpak(jj) /= 0 .and. this%idxbuduzet > 0) then + call this%uze_uzet_term(jj, n1, n2, v) + end if + case ('REJ-INF-TO-MVR') + if (this%iboundpak(jj) /= 0 .and. this%idxbudritm > 0) then + call this%uze_ritm_term(jj, n1, n2, v) + end if + case ('THERMAL-EQUIL') + if (this%iboundpak(jj) /= 0 .and. this%idxbudtheq > 0) then + call this%uze_theq_term(jj, n1, n2, v) + end if + case default + found = .false. + end select + ! + ! -- Return + return + end subroutine uze_bd_obs + + !> @brief Sets the stress period attributes for keyword use. + !< + subroutine uze_set_stressperiod(this, itemno, keyword, found) + ! -- modules + use TimeSeriesManagerModule, only: read_value_or_time_series_adv + ! -- dummy + class(GweUzeType), intent(inout) :: this + integer(I4B), intent(in) :: itemno + character(len=*), intent(in) :: keyword + logical, intent(inout) :: found + ! -- local + character(len=LINELENGTH) :: temp_text + integer(I4B) :: ierr + integer(I4B) :: jj + real(DP), pointer :: bndElem => null() + ! + ! INFILTRATION + ! UZET + ! + found = .true. + select case (keyword) + case ('INFILTRATION') + ierr = this%apt_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(temp_text) + jj = 1 + bndElem => this%tempinfl(itemno) + call read_value_or_time_series_adv(temp_text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'INFILTRATION') + case ('UZET') + ierr = this%apt_check_valid(itemno) + if (ierr /= 0) then + goto 999 + end if + call this%parser%GetString(temp_text) + jj = 1 + bndElem => this%tempuzet(itemno) + call read_value_or_time_series_adv(temp_text, itemno, jj, bndElem, & + this%packName, 'BND', this%tsManager, & + this%iprpak, 'UZET') + case default + ! + ! -- Keyword not recognized so return to caller with found = .false. + found = .false. + end select + ! +999 continue + ! + ! -- Return + return + end subroutine uze_set_stressperiod + +end module GweUzeModule diff --git a/src/Model/TransportModel/tsp1apt1.f90 b/src/Model/TransportModel/tsp1apt1.f90 index 76fc718e033..b63a3b70569 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', 'MWE') + case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE', 'UZE') 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', 'MWE') + case ('LKT', 'SFT', 'MWT', 'UZT', 'LKE', 'SFE', 'MWE', 'UZE') 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 3eca98fdb7e..69ebe64dab6 100644 --- a/src/meson.build +++ b/src/meson.build @@ -73,6 +73,7 @@ modflow_sources = files( 'Model' / 'GroundWaterEnergy' / 'gwe1lke1.f90', 'Model' / 'GroundWaterEnergy' / 'gwe1mwe1.f90', 'Model' / 'GroundWaterEnergy' / 'gwe1sfe1.f90', + 'Model' / 'GroundWaterEnergy' / 'gwe1uze1.f90', 'Model' / 'GroundWaterFlow' / 'gwf3.f90', 'Model' / 'GroundWaterFlow' / 'gwf3api8.f90', 'Model' / 'GroundWaterFlow' / 'gwf3buy8.f90',