diff --git a/autotest/test_snf_disl.py b/autotest/test_snf_disl.py new file mode 100644 index 00000000000..b9121dfdfce --- /dev/null +++ b/autotest/test_snf_disl.py @@ -0,0 +1,186 @@ +""" + +Simple 4 segment network with 5 vertices + + +zero-based diagram below +v1 +o + | + | + | s0 + | + | + | +o------o------o------o +v0 v2 v3 v4 + s1 s2 s3 + +ia ja +0 0 2 +2 1 2 +4 2 0 1 3 +8 3 2 +10 + +""" + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + +ex = ["snf-disl01",] + + +def build_model(idx, dir): + + sim_ws = dir + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=sim_ws, + memory_print_option='all', + ) + + tdis = flopy.mf6.ModflowTdis(sim) + ems = flopy.mf6.ModflowEms(sim) + snf = flopy.mf6.ModflowSnf(sim, modelname=name, save_flows=True) + + vertices = [ + [0, 0., 0., 0.], + [1, 0., 1., 0.], + [2, 1., 0., 0.], + [3, 2., 0., 0.], + [4, 3., 0., 0.], + ] + # icell1d fdc ncvert icvert + cell2d = [ + [0, 0.5, 2, 1, 2], + [1, 0.5, 2, 0, 2], + [2, 0.5, 2, 2, 3], + [3, 0.5, 2, 3, 4], + ] + + nodes = len(cell2d) + nvert = len(vertices) + + disl = flopy.mf6.ModflowSnfdisl( + snf, + nodes=nodes, + nvert=nvert, + segment_length=1000.0, + tosegment=[2, 2, 3, -1], # -1 gives 0 in one-based, which means outflow cell + idomain=1, + vertices=vertices, + cell2d=cell2d, + ) + + # note: for specifying zero-based reach number, put reach number in tuple + fname = f"{name}.mmr.obs.csv" + mmr_obs = { + fname: [ + ("OUTFLOW", "EXT-OUTFLOW", (nodes - 1,)), + ], + "digits": 10, + } + + mmr = flopy.mf6.ModflowSnfmmr( + snf, + observations=mmr_obs, + print_flows=True, + save_flows=True, + iseg_order=list(range(nodes)), + qoutflow0=0.0, + k_coef=0.001, + x_coef=0.2 + ) + + # output control + oc = flopy.mf6.ModflowSnfoc( + snf, + budget_filerecord=f"{name}.bud", + saverecord=[("BUDGET", "ALL"), ], + printrecord=[("BUDGET", "ALL"), ], + ) + + flw_spd = {0: [[0, 1000.0], [1, 500.]]} + flw = flopy.mf6.ModflowSnfflw( + snf, + print_input=True, + print_flows=True, + stress_period_data=flw_spd, + ) + + return sim, None + + +def eval_model(sim): + print("evaluating model...") + + # read the observation output + name = ex[sim.idxsim] + fpth = os.path.join(sim.simpath, f"{name}.mmr.obs.csv") + obsvals = np.genfromtxt(fpth, names=True, delimiter=",") + + # read the binary grid file + fpth = os.path.join(sim.simpath, f"{name}.disl.grb") + grb = flopy.mf6.utils.MfGrdFile(fpth) + ia = grb.ia + ja = grb.ja + assert ia.shape[0] == grb.nodes + 1, "ia in grb file is not correct size" + + # read the budget file + fpth = os.path.join(sim.simpath, f"{name}.bud") + budobj = flopy.utils.binaryfile.CellBudgetFile(fpth) + flowja = budobj.get_data(text="FLOW-JA-FACE") + qstorage = budobj.get_data(text="STORAGE") + qflw = budobj.get_data(text="FLW") + qextoutflow = budobj.get_data(text="EXT-OUTFLOW") + qresidual = np.zeros(grb.nodes) + + # check budget terms + for itime in range(len(flowja)): + print (f"evaluating timestep {itime}") + + fja = flowja[itime].flatten() + for n in range(grb.nodes): + ipos = ia[n] + qresidual[n] = fja[ipos] + assert np.allclose(qresidual, 0.), "residual in flowja diagonal is not zero" + + for n in range(grb.nodes): + qs = qstorage[itime].flatten()[n] + if n + 1 in qflw[itime]["node"]: + idx, = np.where(qflw[itime]["node"] == n + 1) + idx = idx[0] + qf = qflw[itime].flatten()["q"][idx] + else: + qf = 0. + qe = qextoutflow[itime].flatten()[n] + qdiag = fja[ia[n]] + print(f"{n=} {qs=} {qf=} {qe=} {qdiag=}") + for ipos in range(ia[n] + 1, ia[n + 1]): + j = ja[ipos] + q = fja[ipos] + print(f" {ipos=} {j=} {q=}") + + return + + +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + test = TestFramework() + test.build(build_model, idx, ws) + test.run( + TestSimulation( + name=name, exe_dict=targets, exfunc=eval_model, idxsim=idx + ), + ws, + ) diff --git a/autotest/test_snf_muskingum.py b/autotest/test_snf_muskingum.py new file mode 100644 index 00000000000..62f202d7b70 --- /dev/null +++ b/autotest/test_snf_muskingum.py @@ -0,0 +1,187 @@ +""" + +http://uon.sdsu.edu/enghydro/engineering_hydrology_09.php#muskingum +Ponce -- Muskingum routing example +single reach +""" + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + + +# data from ponce Table 9-1 +# time(d) coi2(m4/s) c1i1(m3/s) c2i1(m3/s) outflow(m3/s) +ponce_data = """0 352.0 0 0 0 352.0 +1 587.0 76.6 107.1 199.0 382.7 +2 1353.0 176.5 178.6 216.3 571.4 +3 2725.0 355.4 411.8 323.0 1090.2 +4 4408.5 575.0 829.4 616.2 2020.6 +5 5987.0 780.9 1341.7 1142.1 3264.7 +6 6704.0 874.4 1822.1 1845.3 4541.8 +7 6951.0 906.7 2040.3 2567.1 5514.1 +8 6839.0 892.0 2115.5 3116.7 6124.2 +9 6207.0 809.6 2081.5 3461.5 6352.6 +10 5346.0 697.3 1889.1 3590.6 6177.0 +11 4560.0 594.8 1627.0 3491.4 5713.2 +12 3861.5 503.7 1387.8 3229.2 5120.7 +13 3007.0 392.2 1175.2 2894.3 4461.7 +14 2357.5 307.5 915.2 2521.8 3744.5 +15 1779.0 232.0 717.5 2116.5 3066.0 +16 1405.0 183.3 541.4 1733.0 2457.7 +17 1123.0 146.5 427.6 1389.1 1963.2 +18 952.5 124.2 341.8 1109.6 1575.6 +19 730.0 95.2 289.9 890.6 1275.7 +20 605.0 78.9 222.2 721.0 1022.1 +21 514.0 67.1 184.1 577.7 828.9 +22 422.0 55.1 156.4 468.5 680.0 +23 352.0 45.9 128.4 384.4 558.7 +24 352.0 45.9 107.1 315.8 468.8 +25 352.0 45.9 107.1 265.0 418.0""" + + +def get_ponce_data(): + qinflow = [] + qoutflow = [] + time_days = [] + for line in ponce_data.split("\n"): + itime, inflow, c0, c1, c2, outflow = line.strip().split(" ") + time_days.append(float(itime)) + qinflow.append(float(inflow)) + qoutflow.append(float(outflow)) + return time_days, qinflow, qoutflow + + +ex = ["ponce01",] + + +def build_model(idx, dir): + + sim_ws = dir + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=sim_ws, + memory_print_option='all', + ) + nper = 25 + tdis_rc = [(1., 1, 1.0) for ispd in range(nper)] + tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_rc) + ems = flopy.mf6.ModflowEms(sim) + snf = flopy.mf6.ModflowSnf(sim, modelname=name, save_flows=True) + + vertices = None + cell2d = None + nvert = None + + nodes = 1 + channel_length = 500. # miles + channel_length = channel_length * 5280. / 3.2808 # meters + segment_length = channel_length * 5280. / nodes + + k_coef = 2. # days + x_coef = 0.1 # dimensionless + + time_days, qinflow, qoutflow = get_ponce_data() + + disl = flopy.mf6.ModflowSnfdisl( + snf, + nodes=nodes, + nvert=nvert, + segment_length=segment_length, + tosegment=[(-1,)], # (-1,) gives 0 in one-based, which means outflow cell + idomain=1, + vertices=vertices, + cell2d=cell2d, + ) + + + # note: for specifying reach number, use fortran indexing! + fname = f"{name}.mmr.obs.csv" + mmr_obs = { + fname: [ + ("OUTFLOW1", "EXT-OUTFLOW", 1), + ], + "digits": 10, + } + + mmr = flopy.mf6.ModflowSnfmmr( + snf, + print_flows=True, + observations=mmr_obs, + iseg_order=list(range(nodes)), + qoutflow0=qinflow[0], + k_coef=k_coef, + x_coef=x_coef + ) + + # output control + oc = flopy.mf6.ModflowSnfoc( + snf, + budget_filerecord=f"{name}.bud", + saverecord=[("BUDGET", "ALL"), ], + printrecord=[("BUDGET", "ALL"), ], + ) + + inflow = qinflow[1:] + flw_spd = {ispd: [[0, inflow[ispd]]] for ispd in range(nper)} + flw = flopy.mf6.ModflowSnfflw( + snf, + print_input=True, + print_flows=True, + stress_period_data=flw_spd, + ) + + return sim, None + + +def eval_model(sim): + print("evaluating model...") + + # get back the ponce data for comparison + time_days, qinflow, qoutflow = get_ponce_data() + + # read the observation output + name = ex[sim.idxsim] + fpth = os.path.join(sim.simpath, f"{name}.mmr.obs.csv") + obsvals = np.genfromtxt(fpth, names=True, delimiter=",") + + # compare output with known result + time_days = time_days[1:] + qoutflow = -np.array(qoutflow[1:]) + atol = 0.06 + success = np.allclose(qoutflow, obsvals["OUTFLOW1"], atol=atol) + if not success: + for i, t in enumerate(time_days): + qa = qoutflow[i] + qs = obsvals["OUTFLOW1"][i] + d = qa - qs + if i == 0: + maxdiff = d + else: + if abs(d) > maxdiff: + maxdiff = abs(d) + print(t, qa, qs, d) + print(f"maximum difference is {maxdiff}") + assert success + + return + + +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + test = TestFramework() + test.build(build_model, idx, ws) + test.run( + TestSimulation( + name=name, exe_dict=targets, exfunc=eval_model, idxsim=idx + ), + ws, + ) diff --git a/autotest/test_snf_thomas.py b/autotest/test_snf_thomas.py new file mode 100644 index 00000000000..833fdd4935e --- /dev/null +++ b/autotest/test_snf_thomas.py @@ -0,0 +1,263 @@ +""" + +http://uon.sdsu.edu/the_thomas_problem_with_online_computation.html +Thomas routine example +500 miles +dx = 25 miles +sinusoidal inflow + +""" + +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + + +ponce_data = """0 0 50 50 +1 3 51.441 50 +2 6 55.709 50 +3 9 62.64 50 +4 12 71.967 50 +5 15 83.332 50 +6 18 96.299 50 +7 21 110.368 50 +8 24 125 50 +9 27 139.632 50 +10 30 153.701 50 +11 33 166.668 50 +12 36 178.033 50 +13 39 187.36 50 +14 42 194.291 50 +15 45 198.559 50.002 +16 48 200 50.006 +17 51 198.559 50.018 +18 54 194.291 50.049 +19 57 187.36 50.121 +20 60 178.033 50.275 +21 63 166.668 50.575 +22 66 153.701 51.119 +23 69 139.632 52.034 +24 72 125 53.476 +25 75 110.368 55.616 +26 78 96.299 58.627 +27 81 83.332 62.656 +28 84 71.967 67.811 +29 87 62.64 74.142 +30 90 55.709 81.625 +31 93 51.441 90.165 +32 96 50 99.589 +33 99 50 109.662 +34 102 50 120.091 +35 105 50 130.549 +36 108 50 140.685 +37 111 50 150.146 +38 114 50 158.595 +39 117 50 165.722 +40 120 50 171.266 +41 123 50 175.021 +42 126 50 176.847 +43 129 50 176.676 +44 132 50 174.517 +45 135 50 170.453 +46 138 50 164.641 +47 141 50 157.305 +48 144 50 148.724 +49 147 50 139.223 +50 150 50 129.156 +51 153 50 118.886 +52 156 50 108.765 +53 159 50 99.11 +54 162 50 90.185 +55 165 50 82.184 +56 168 50 75.223 +57 171 50 69.342 +58 174 50 64.515 +59 177 50 60.662 +60 180 50 57.669 +61 183 50 55.403 +62 186 50 53.732 +63 189 50 52.528 +64 192 50 51.68 +65 195 50 51.097 +66 198 50 50.704 +67 201 50 50.444 +68 204 50 50.275 +69 207 50 50.168 +70 210 50 50.101 +71 213 50 50.06 +72 216 50 50.035 +73 219 50 50.02 +74 222 50 50.012 +75 225 50 50.007 +76 228 50 50.004 +77 231 50 50.002 +78 234 50 50.001 +79 237 50 50.001 +80 240 50 50""" + + +def get_ponce_data(): + qinflow = [] + qoutflow = [] + time_days = [] + for line in ponce_data.split("\n"): + itime, elapsed_time, inflow, outflow = line.strip().split(" ") + time_days.append(float(elapsed_time)) + qinflow.append(float(inflow)) + qoutflow.append(float(outflow)) + return time_days, qinflow, qoutflow + + +ex = ["thomas01",] + + +def build_model(idx, dir): + + sim_ws = dir + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=sim_ws, + memory_print_option='all', + ) + delta_t = 60 * 60 * 3 # seconds = 3 hours + nper = 80 + tdis_rc = [(delta_t, 1, 1.0) for ispd in range(nper)] + tdis = flopy.mf6.ModflowTdis(sim, nper=nper, perioddata=tdis_rc, time_units="seconds") + ems = flopy.mf6.ModflowEms(sim) + snf = flopy.mf6.ModflowSnf(sim, modelname=name, save_flows=True) + + vertices = None + cell2d = None + nvert = None + + nodes = 20 # can run 200,000 reaches in less than a second and 2 million reaches in 4 seconds + channel_length = 500 * 5280. # meters + reach_length = channel_length / nodes + + qp = 200 # cfs = peak discharge + qb = 50 # cfs = base discharge + qa = (qp + qb) / 2 + + beta = 5. / 3. + rating_coefficient = 0.688 + depth_a = 2. + # solve for reference depth from rating equation + # reference depth is the depth when flow is reference flow qa + depth_a = (qa / rating_coefficient) ** (1 / beta) + v_mean = qa / depth_a # mean velocity at peak discharge + wave_celerity = beta * v_mean # 4.0 # m/s calculated as beta * v + q0 = qa # flow per unit width based on average discharge + S0 = 1. / 5280. # channel bottom slope, 1 ft/mi converted to ft/ft + + print(f"{depth_a=} ft") + print(f"{v_mean=} ft/s") + print(f"{wave_celerity=} ft/s") + + k_coef = reach_length / wave_celerity # seconds + x_coef = 0.5 * (1 - q0 / S0 / wave_celerity / reach_length) # dimensionless + + time_days, qinflow, qoutflow = get_ponce_data() + + tosegment = [(irch,) for irch in range(1, nodes)] + [(-1,)] + disl = flopy.mf6.ModflowSnfdisl( + snf, + length_units="feet", + nodes=nodes, + nvert=nvert, + segment_length=reach_length, + tosegment=tosegment, # (-1,) gives 0 in one-based, which means outflow cell + idomain=1, + vertices=vertices, + cell2d=cell2d, + ) + + # note: for specifying reach number, use fortran indexing! + fname = f"{name}.mmr.obs.csv" + mmr_obs = { + fname: [ + ("OUTFLOW", "EXT-OUTFLOW", (nodes - 1,)), + ], + "digits": 10, + } + + mmr = flopy.mf6.ModflowSnfmmr( + snf, + print_flows=False, + observations=mmr_obs, + iseg_order=list(range(nodes)), + qoutflow0=qinflow[0], + k_coef=k_coef, + x_coef=x_coef + ) + + # output control + oc = flopy.mf6.ModflowSnfoc( + snf, + budget_filerecord=f"{name}.bud", + saverecord=[("BUDGET", "ALL"), ], + printrecord=[("BUDGET", "ALL"), ], + ) + + inflow = qinflow[1:] + flw_spd = {ispd: [[0, inflow[ispd]]] for ispd in range(nper)} + flw = flopy.mf6.ModflowSnfflw( + snf, + print_input=True, + print_flows=True, + stress_period_data=flw_spd, + ) + + return sim, None + + +def eval_model(sim): + print("evaluating model...") + + # get back the ponce data for comparison + time_days, qinflow, qoutflow = get_ponce_data() + + # read the observation output + name = ex[sim.idxsim] + fpth = os.path.join(sim.simpath, f"{name}.mmr.obs.csv") + obsvals = np.genfromtxt(fpth, names=True, delimiter=",") + + # compare output with known result + time_days = time_days[1:] + qoutflow = -np.array(qoutflow[1:]) + atol = 0.001 + success = np.allclose(qoutflow, obsvals["OUTFLOW"], atol=atol) + if not success: + for i, t in enumerate(time_days): + qa = qoutflow[i] + qs = obsvals["OUTFLOW"][i] + d = qa - qs + if i == 0: + maxdiff = d + else: + if abs(d) > maxdiff: + maxdiff = abs(d) + print(t, qa, qs, d) + print(f"maximum difference is {maxdiff}") + assert success + + return + + +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + test = TestFramework() + test.build(build_model, idx, ws) + test.run( + TestSimulation( + name=name, exe_dict=targets, exfunc=eval_model, idxsim=idx + ), + ws, + ) diff --git a/doc/mf6io/mf6ivar/dfn/sim-nam.dfn b/doc/mf6io/mf6ivar/dfn/sim-nam.dfn index ac9be58e7bb..fda8e38fbe4 100644 --- a/doc/mf6io/mf6ivar/dfn/sim-nam.dfn +++ b/doc/mf6io/mf6ivar/dfn/sim-nam.dfn @@ -163,7 +163,7 @@ description is the list of solution types and models in the solution. block solutiongroup name slntype type string -valid ims6 +valid ims6 ems6 in_record true tagged false reader urword diff --git a/doc/mf6io/mf6ivar/dfn/sln-ems.dfn b/doc/mf6io/mf6ivar/dfn/sln-ems.dfn new file mode 100644 index 00000000000..03ff8e88694 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/sln-ems.dfn @@ -0,0 +1 @@ +# flopy solution_package ems * diff --git a/doc/mf6io/mf6ivar/dfn/snf-disl.dfn b/doc/mf6io/mf6ivar/dfn/snf-disl.dfn new file mode 100644 index 00000000000..9e365def4f5 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/snf-disl.dfn @@ -0,0 +1,218 @@ +# --------------------- snf disl options --------------------- + +block options +name length_units +type string +reader urword +optional true +longname model length units +description is the length units used for this model. Values can be ``FEET'', ``METERS'', or ``CENTIMETERS''. If not specified, the default is ``UNKNOWN''. + +block options +name length_convert +type double precision +reader urword +optional true +longname length conversion factor +description conversion factor that is used in converting constants used by MODFLOW into model length units. All constants are by default in units of "meters". Constants that use length_conversion include GRAVITY and STORAGE COEFFICIENT values. + +block options +name time_convert +type double precision +reader urword +optional true +longname time conversion factor +description conversion factor that is used in converting constants used by MODFLOW into model time units. All constants are by default in units of "seconds". TIME_CONVERSION is used to convert the GRAVITY constant from seconds to the model's units. + +block options +name nogrb +type keyword +reader urword +optional true +longname do not write binary grid file +description keyword to deactivate writing of the binary grid file. + +block options +name xorigin +type double precision +reader urword +optional true +longname x-position origin of the model grid coordinate system +description x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +block options +name yorigin +type double precision +reader urword +optional true +longname y-position origin of the model grid coordinate system +description y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +block options +name angrot +type double precision +reader urword +optional true +longname rotation angle +description counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +# --------------------- snf disl dimensions --------------------- + +block dimensions +name nodes +type integer +reader urword +optional false +longname number of linear features +description is the number of linear cells. + +block dimensions +name nvert +type integer +reader urword +optional true +longname number of columns +description is the total number of (x, y, z) vertex pairs used to characterize the model grid. + +# --------------------- snf disl griddata --------------------- + +block griddata +name segment_length +type double precision +shape (nodes) +valid +reader readarray +layered false +optional +longname segment length +description length for each segment + +block griddata +name tosegment +type integer +shape (nodes) +valid +reader readarray +layered false +optional +longname index of the downstream segment +description index of downstream segment to which the segment streamflow flows, for segments that do not flow to another segment enter 0 +numeric_index true + +block griddata +name idomain +type integer +shape (nodes) +reader readarray +layered false +optional true +longname idomain existence array +description is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. + + +# --------------------- snf disl vertices --------------------- + +block vertices +name vertices +type recarray iv xv yv zv +shape (nvert) +reader urword +optional false +longname vertices data +description + +block vertices +name iv +type integer +in_record true +tagged false +reader urword +optional false +longname vertex number +description is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. +numeric_index true + +block vertices +name xv +type double precision +in_record true +tagged false +reader urword +optional false +longname x-coordinate for vertex +description is the x-coordinate for the vertex. + +block vertices +name yv +type double precision +in_record true +tagged false +reader urword +optional false +longname y-coordinate for vertex +description is the y-coordinate for the vertex. + +block vertices +name zv +type double precision +in_record true +tagged false +reader urword +optional false +longname z-coordinate for vertex +description is the z-coordinate for the vertex. + + +# --------------------- snf disl cell2d --------------------- + +block cell2d +name cell2d +type recarray icell2d fdc ncvert icvert +shape (nodes) +reader urword +optional false +longname cell2d data +description + +block cell2d +name icell2d +type integer +in_record true +tagged false +reader urword +optional false +longname cell2d number +description is the cell2d number. Records in the cell2d block must be listed in consecutive order from the first to the last. +numeric_index true + +block cell2d +name fdc +type double precision +in_record true +tagged false +reader urword +optional false +longname fractional distance to the cell center +description is the fractional distance to the cell center. FDC is relative to the first vertex in the ICVERT array. In most cases FDC should be 0.5, which would place the center of the line segment that defines the cell. If the value of FDC is 1, the cell center would located at the last vertex. FDC values of 0 and 1 can be used to place the node at either end of the cell which can be useful for cells with boundary conditions. + +block cell2d +name ncvert +type integer +in_record true +tagged false +reader urword +optional false +longname number of cell vertices +description is the number of vertices required to define the cell. There may be a different number of vertices for each cell. + +block cell2d +name icvert +type integer +shape (ncvert) +in_record true +tagged false +reader urword +optional false +longname number of cell vertices +description is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in the order that defines the line representing the cell. Cells that are connected must share vertices. The bottom elevation of the cell is calculated using the ZV of the first and last vertex point and FDC. +numeric_index true \ No newline at end of file diff --git a/doc/mf6io/mf6ivar/dfn/snf-flw.dfn b/doc/mf6io/mf6ivar/dfn/snf-flw.dfn new file mode 100644 index 00000000000..7afba843299 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/snf-flw.dfn @@ -0,0 +1,206 @@ +# --------------------- snf flw options --------------------- +# flopy multi-package + +block options +name auxiliary +type string +shape (naux) +reader urword +optional true +longname keyword to specify aux variables +description REPLACE auxnames {'{#1}': 'Stream Network Flow'} + +block options +name auxmultname +type string +shape +reader urword +optional true +longname name of auxiliary variable for multiplier +description REPLACE auxmultname {'{#1}': 'flow rate'} + +block options +name boundnames +type keyword +shape +reader urword +optional true +longname +description REPLACE boundnames {'{#1}': 'inflow'} + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'inflow'} + +block options +name print_flows +type keyword +reader urword +optional true +longname print calculated flows to listing file +description REPLACE print_flows {'{#1}': 'inflow'} + +block options +name save_flows +type keyword +reader urword +optional true +longname save well flows to budget file +description REPLACE save_flows {'{#1}': 'inflow'} + +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}': 'inflow'} + +# --------------------- snf flw dimensions --------------------- + +block dimensions +name maxbound +type integer +reader urword +optional false +longname maximum number of inflow +description REPLACE maxbound {'{#1}': 'inflow'} + + +# --------------------- snf flw 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 stress_period_data +type recarray cellid q aux boundname +shape (maxbound) +reader urword +longname +description + +block period +name cellid +type integer +shape (ncelldim) +tagged false +in_record true +reader urword +longname cell identifier +description REPLACE cellid {} + +block period +name q +type double precision +shape +tagged false +in_record true +reader urword +time_series true +longname well rate +description is the volumetric inflow rate. A positive value indicates inflow to the stream. Negative values are not allows. 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 aux +type double precision +in_record true +tagged false +shape (naux) +reader urword +optional true +time_series true +longname auxiliary variables +description REPLACE aux {'{#1}': 'inflow'} + +block period +name boundname +type string +shape +tagged false +in_record true +reader urword +optional true +longname inflow name +description REPLACE boundname {'{#1}': 'inflow'} diff --git a/doc/mf6io/mf6ivar/dfn/snf-mmr.dfn b/doc/mf6io/mf6ivar/dfn/snf-mmr.dfn new file mode 100644 index 00000000000..4ee6dffb412 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/snf-mmr.dfn @@ -0,0 +1,108 @@ +# --------------------- snf mmr options --------------------- + +block options +name save_flows +type keyword +reader urword +optional true +longname keyword to save NPF flows +description keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. +mf6internal ipakcb + +block options +name print_flows +type keyword +reader urword +optional true +longname keyword to print MMR flows to listing file +description keyword to indicate that calculated flows between cells 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. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the MMR Package for models that have a small number of cells. +mf6internal iprflow + +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 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 obs6_filename +type string +in_record true +tagged false +reader urword +optional false +longname obs6 input filename +description REPLACE obs6_filename {'{#1}': 'MMR'} + +# --------------------- snf mmr griddata --------------------- + +block griddata +name iseg_order +type integer +shape (nodes) +valid +reader readarray +layered false +optional +longname segment calculation order +description segment calculation order +numeric_index true + +block griddata +name qoutflow0 +type double precision +shape (nodes) +valid +reader readarray +layered false +optional +longname initial flows +description initial outflow from the reach at time zero + +block griddata +name k_coef +type double precision +shape (nodes) +valid +reader readarray +layered false +optional +longname manning K coefficient +description manning K coefficient + +block griddata +name x_coef +type double precision +shape (nodes) +valid +reader readarray +layered false +optional +longname routing weighting factor +description The amount of attenuation of the flow wave, called the Muskingum routing weighting factor; enter 0.0 for reservoirs, diversions, and segment(s) flowing out of the basin diff --git a/doc/mf6io/mf6ivar/dfn/snf-nam.dfn b/doc/mf6io/mf6ivar/dfn/snf-nam.dfn new file mode 100644 index 00000000000..7a26000dcfc --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/snf-nam.dfn @@ -0,0 +1,73 @@ +# --------------------- snf nam options --------------------- + +block options +name list +type string +reader urword +optional true +longname name of listing file +description is name of the listing file to create for this SNF model. If not specified, then the name of the list file will be the basename of the SNF model name file and the '.lst' extension. For example, if the SNF name file is called ``my.model.nam'' then the list file will be called ``my.model.lst''. + +block options +name print_input +type keyword +reader urword +optional true +longname print input to listing file +description REPLACE print_input {'{#1}': 'all model stress package'} + +block options +name print_flows +type keyword +reader urword +optional true +longname print calculated flows to listing file +description REPLACE print_flows {'{#1}': 'all model package'} + +block options +name save_flows +type keyword +reader urword +optional true +longname save flows for all packages to budget file +description REPLACE save_flows {'{#1}': 'all model package'} + +# --------------------- snf nam packages --------------------- + +block packages +name packages +type recarray ftype fname pname +reader urword +optional false +longname package list +description + +block packages +name ftype +in_record true +type string +tagged false +reader urword +longname package type +description is the file type, which must be one of the following character values shown in table~\ref{table:ftype}. Ftype may be entered in any combination of uppercase and lowercase. + +block packages +name fname +in_record true +type string +preserve_case true +tagged false +reader urword +longname file name +description is the name of the file containing the package input. The path to the file should be included if the file is not located in the folder where the program was run. + +block packages +name pname +in_record true +type string +tagged false +reader urword +optional true +longname user name for package +description is the user-defined name for the package. PNAME is restricted to 16 characters. No spaces are allowed in PNAME. PNAME character values are read and stored by the program for stress packages only. These names may be useful for labeling purposes when multiple stress packages of the same type are located within a single SNF Model. If PNAME is specified for a stress package, then PNAME will be used in the flow budget table in the listing file; it will also be used for the text entry in the cell-by-cell budget file. PNAME is case insensitive and is stored in all upper case letters. + diff --git a/doc/mf6io/mf6ivar/dfn/snf-oc.dfn b/doc/mf6io/mf6ivar/dfn/snf-oc.dfn new file mode 100644 index 00000000000..ae425b24453 --- /dev/null +++ b/doc/mf6io/mf6ivar/dfn/snf-oc.dfn @@ -0,0 +1,204 @@ +# --------------------- snf oc options --------------------- + +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 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. + +# --------------------- snf oc 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 saverecord +type record save rtype ocsetting +shape +reader urword +tagged false +optional true +longname +description + +block period +name save +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname keyword to save +description keyword to indicate that information will be saved this stress period. + +block period +name printrecord +type record print rtype ocsetting +shape +reader urword +tagged false +optional true +longname +description + +block period +name print +type keyword +shape +in_record true +reader urword +tagged true +optional false +longname keyword to save +description keyword to indicate that information will be printed this stress period. + +block period +name rtype +type string +shape +in_record true +reader urword +tagged false +optional false +longname record type +description type of information to save or print. Can be BUDGET. + +block period +name ocsetting +type keystring all first last frequency steps +shape +tagged false +in_record true +reader urword +longname +description specifies the steps for which the data will be saved. + +block period +name all +type keyword +shape +in_record true +reader urword +longname +description keyword to indicate save for all time steps in period. + +block period +name first +type keyword +shape +in_record true +reader urword +longname +description keyword to indicate save for first step in period. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. + +block period +name last +type keyword +shape +in_record true +reader urword +longname +description keyword to indicate save for last step in period. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. + +block period +name frequency +type integer +shape +tagged true +in_record true +reader urword +longname +description save at the specified time step frequency. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. + +block period +name steps +type integer +shape ($0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect. When THICKSTRT is in effect, a negative value of icelltype indicates that saturated thickness will be computed as STRT-BOT and held constant. | +| GWF | NPF | GRIDDATA | ICELLTYPE | INTEGER (NODES) | flag for each cell that specifies how saturated thickness is treated. 0 means saturated thickness is held constant; $>$0 means saturated thickness varies with computed head when head is below the cell top; $<$0 means saturated thickness varies with computed head unless the THICKSTRT option is in effect. When THICKSTRT is in effect, a negative value for ICELLTYPE indicates that the saturated thickness value used in conductance calculations in the NPF Package will be computed as STRT-BOT and held constant. If the THICKSTRT option is not in effect, then negative values provided by the user for ICELLTYPE are automatically reassigned by the program to a value of one. | | GWF | NPF | GRIDDATA | K | DOUBLE PRECISION (NODES) | is the hydraulic conductivity. For the common case in which the user would like to specify the horizontal hydraulic conductivity and the vertical hydraulic conductivity, then K should be assigned as the horizontal hydraulic conductivity, K33 should be assigned as the vertical hydraulic conductivity, and K22 and the three rotation angles should not be specified. When more sophisticated anisotropy is required, then K corresponds to the K11 hydraulic conductivity axis. All included cells (IDOMAIN $>$ 0) must have a K value greater than zero. | | GWF | NPF | GRIDDATA | K22 | DOUBLE PRECISION (NODES) | is the hydraulic conductivity of the second ellipsoid axis (or the ratio of K22/K if the K22OVERK option is specified); for an unrotated case this is the hydraulic conductivity in the y direction. If K22 is not included in the GRIDDATA block, then K22 is set equal to K. For a regular MODFLOW grid (DIS Package is used) in which no rotation angles are specified, K22 is the hydraulic conductivity along columns in the y direction. For an unstructured DISU grid, the user must assign principal x and y axes and provide the angle for each cell face relative to the assigned x direction. All included cells (IDOMAIN $>$ 0) must have a K22 value greater than zero. | | GWF | NPF | GRIDDATA | K33 | DOUBLE PRECISION (NODES) | is the hydraulic conductivity of the third ellipsoid axis (or the ratio of K33/K if the K33OVERK option is specified); for an unrotated case, this is the vertical hydraulic conductivity. When anisotropy is applied, K33 corresponds to the K33 tensor component. All included cells (IDOMAIN $>$ 0) must have a K33 value greater than zero. | @@ -1151,6 +1151,74 @@ | GWT | API | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the api boundary 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 api boundary package. | | GWT | API | OPTIONS | MOVER | KEYWORD | keyword to indicate that this instance of the api boundary Package can be used with the Water Mover (MVR) Package. When the MOVER option is specified, additional memory is allocated within the package to store the available, provided, and received water. | | GWT | API | DIMENSIONS | MAXBOUND | INTEGER | integer value specifying the maximum number of api boundary cells that will be specified for use during any stress period. | +| SNF | NAM | OPTIONS | LIST | STRING | is name of the listing file to create for this SNF model. If not specified, then the name of the list file will be the basename of the SNF model name file and the '.lst' extension. For example, if the SNF name file is called ``my.model.nam'' then the list file will be called ``my.model.lst''. | +| SNF | NAM | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of all model stress package information will be written to the listing file immediately after it is read. | +| SNF | NAM | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of all model package flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. | +| SNF | NAM | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that all model package flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | +| SNF | NAM | PACKAGES | FTYPE | STRING | is the file type, which must be one of the following character values shown in table~\ref{table:ftype}. Ftype may be entered in any combination of uppercase and lowercase. | +| SNF | NAM | PACKAGES | FNAME | STRING | is the name of the file containing the package input. The path to the file should be included if the file is not located in the folder where the program was run. | +| SNF | NAM | PACKAGES | PNAME | STRING | is the user-defined name for the package. PNAME is restricted to 16 characters. No spaces are allowed in PNAME. PNAME character values are read and stored by the program for stress packages only. These names may be useful for labeling purposes when multiple stress packages of the same type are located within a single SNF Model. If PNAME is specified for a stress package, then PNAME will be used in the flow budget table in the listing file; it will also be used for the text entry in the cell-by-cell budget file. PNAME is case insensitive and is stored in all upper case letters. | +| SNF | DISL | OPTIONS | LENGTH_UNITS | STRING | is the length units used for this model. Values can be ``FEET'', ``METERS'', or ``CENTIMETERS''. If not specified, the default is ``UNKNOWN''. | +| SNF | DISL | OPTIONS | LENGTH_CONVERT | DOUBLE PRECISION | conversion factor that is used in converting constants used by MODFLOW into model length units. All constants are by default in units of "meters". Constants that use length_conversion include GRAVITY and STORAGE COEFFICIENT values. | +| SNF | DISL | OPTIONS | TIME_CONVERT | DOUBLE PRECISION | conversion factor that is used in converting constants used by MODFLOW into model time units. All constants are by default in units of "seconds". TIME_CONVERSION is used to convert the GRAVITY constant from seconds to the model's units. | +| SNF | DISL | OPTIONS | NOGRB | KEYWORD | keyword to deactivate writing of the binary grid file. | +| SNF | DISL | OPTIONS | XORIGIN | DOUBLE PRECISION | x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SNF | DISL | OPTIONS | YORIGIN | DOUBLE PRECISION | y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SNF | DISL | OPTIONS | ANGROT | DOUBLE PRECISION | counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. | +| SNF | DISL | DIMENSIONS | NODES | INTEGER | is the number of linear cells. | +| SNF | DISL | DIMENSIONS | NVERT | INTEGER | is the total number of (x, y, z) vertex pairs used to characterize the model grid. | +| SNF | DISL | GRIDDATA | SEGMENT_LENGTH | DOUBLE PRECISION (NODES) | length for each segment | +| SNF | DISL | GRIDDATA | TOSEGMENT | INTEGER (NODES) | index of downstream segment to which the segment streamflow flows, for segments that do not flow to another segment enter 0 | +| SNF | DISL | GRIDDATA | IDOMAIN | INTEGER (NODES) | is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. | +| SNF | DISL | VERTICES | IV | INTEGER | is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. | +| SNF | DISL | VERTICES | XV | DOUBLE PRECISION | is the x-coordinate for the vertex. | +| SNF | DISL | VERTICES | YV | DOUBLE PRECISION | is the y-coordinate for the vertex. | +| SNF | DISL | VERTICES | ZV | DOUBLE PRECISION | is the z-coordinate for the vertex. | +| SNF | DISL | CELL2D | ICELL2D | INTEGER | is the cell2d number. Records in the cell2d block must be listed in consecutive order from the first to the last. | +| SNF | DISL | CELL2D | FDC | DOUBLE PRECISION | is the fractional distance to the cell center. FDC is relative to the first vertex in the ICVERT array. In most cases FDC should be 0.5, which would place the center of the line segment that defines the cell. If the value of FDC is 1, the cell center would located at the last vertex. FDC values of 0 and 1 can be used to place the node at either end of the cell which can be useful for cells with boundary conditions. | +| SNF | DISL | CELL2D | NCVERT | INTEGER | is the number of vertices required to define the cell. There may be a different number of vertices for each cell. | +| SNF | DISL | CELL2D | ICVERT | INTEGER (NCVERT) | is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in the order that defines the line representing the cell. Cells that are connected must share vertices. The bottom elevation of the cell is calculated using the ZV of the first and last vertex point and FDC. | +| SNF | MMR | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. | +| SNF | MMR | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that calculated flows between cells 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. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the MMR Package for models that have a small number of cells. | +| SNF | MMR | OPTIONS | OBS6 | KEYWORD | keyword to specify that record corresponds to an observations file. | +| SNF | MMR | OPTIONS | FILEIN | KEYWORD | keyword to specify that an input filename is expected next. | +| SNF | MMR | OPTIONS | OBS6_FILENAME | STRING | name of input file to define observations for the MMR 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 MMR package. | +| SNF | MMR | GRIDDATA | ISEG_ORDER | INTEGER (NODES) | segment calculation order | +| SNF | MMR | GRIDDATA | QOUTFLOW0 | DOUBLE PRECISION (NODES) | initial outflow from the reach at time zero | +| SNF | MMR | GRIDDATA | K_COEF | DOUBLE PRECISION (NODES) | manning K coefficient | +| SNF | MMR | GRIDDATA | X_COEF | DOUBLE PRECISION (NODES) | The amount of attenuation of the flow wave, called the Muskingum routing weighting factor; enter 0.0 for reservoirs, diversions, and segment(s) flowing out of the basin | +| SNF | OC | OPTIONS | BUDGET | KEYWORD | keyword to specify that record corresponds to the budget. | +| SNF | OC | OPTIONS | FILEOUT | KEYWORD | keyword to specify that an output filename is expected next. | +| SNF | OC | OPTIONS | BUDGETFILE | STRING | name of the output file to write budget information. | +| SNF | OC | OPTIONS | BUDGETCSV | KEYWORD | keyword to specify that record corresponds to the budget CSV. | +| SNF | OC | 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. | +| SNF | OC | 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. | +| SNF | OC | PERIOD | SAVE | KEYWORD | keyword to indicate that information will be saved this stress period. | +| SNF | OC | PERIOD | PRINT | KEYWORD | keyword to indicate that information will be printed this stress period. | +| SNF | OC | PERIOD | RTYPE | STRING | type of information to save or print. Can be BUDGET. | +| SNF | OC | PERIOD | OCSETTING | KEYSTRING | specifies the steps for which the data will be saved. | +| SNF | OC | PERIOD | ALL | KEYWORD | keyword to indicate save for all time steps in period. | +| SNF | OC | PERIOD | FIRST | KEYWORD | keyword to indicate save for first step in period. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. | +| SNF | OC | PERIOD | LAST | KEYWORD | keyword to indicate save for last step in period. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. | +| SNF | OC | PERIOD | FREQUENCY | INTEGER | save at the specified time step frequency. This keyword may be used in conjunction with other keywords to print or save results for multiple time steps. | +| SNF | OC | PERIOD | STEPS | INTEGER ( + + ... +END CELL1D diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-cell2d.dat b/doc/mf6io/mf6ivar/tex/snf-disl-cell2d.dat new file mode 100644 index 00000000000..f8ad1f65f50 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-cell2d.dat @@ -0,0 +1,5 @@ +BEGIN CELL2D + + + ... +END CELL2D diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-desc.tex b/doc/mf6io/mf6ivar/tex/snf-disl-desc.tex new file mode 100644 index 00000000000..bfbd1a034c0 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-desc.tex @@ -0,0 +1,63 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\item \texttt{length\_units}---is the length units used for this model. Values can be ``FEET'', ``METERS'', or ``CENTIMETERS''. If not specified, the default is ``UNKNOWN''. + +\item \texttt{length\_convert}---conversion factor that is used in converting constants used by MODFLOW into model length units. All constants are by default in units of "meters". Constants that use length_conversion include GRAVITY and STORAGE COEFFICIENT values. + +\item \texttt{time\_convert}---conversion factor that is used in converting constants used by MODFLOW into model time units. All constants are by default in units of "seconds". TIME_CONVERSION is used to convert the GRAVITY constant from seconds to the model's units. + +\item \texttt{NOGRB}---keyword to deactivate writing of the binary grid file. + +\item \texttt{xorigin}---x-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. A default value of zero is assigned if not specified. The value for XORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +\item \texttt{yorigin}---y-position of the origin used for model grid vertices. This value should be provided in a real-world coordinate system. If not specified, then a default value equal to zero is used. The value for YORIGIN does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +\item \texttt{angrot}---counter-clockwise rotation angle (in degrees) of the model grid coordinate system relative to a real-world coordinate system. If not specified, then a default value of 0.0 is assigned. The value for ANGROT does not affect the model simulation, but it is written to the binary grid file so that postprocessors can locate the grid in space. + +\end{description} +\item \textbf{Block: DIMENSIONS} + +\begin{description} +\item \texttt{nodes}---is the number of linear cells. + +\item \texttt{nvert}---is the total number of (x, y, z) vertex pairs used to characterize the model grid. + +\end{description} +\item \textbf{Block: GRIDDATA} + +\begin{description} +\item \texttt{segment\_length}---length for each segment + +\item \texttt{tosegment}---index of downstream segment to which the segment streamflow flows, for segments that do not flow to another segment enter 0 + +\item \texttt{idomain}---is an optional array that characterizes the existence status of a cell. If the IDOMAIN array is not specified, then all model cells exist within the solution. If the IDOMAIN value for a cell is 0, the cell does not exist in the simulation. Input and output values will be read and written for the cell, but internal to the program, the cell is excluded from the solution. If the IDOMAIN value for a cell is 1, the cell exists in the simulation. + +\end{description} +\item \textbf{Block: VERTICES} + +\begin{description} +\item \texttt{iv}---is the vertex number. Records in the VERTICES block must be listed in consecutive order from 1 to NVERT. + +\item \texttt{xv}---is the x-coordinate for the vertex. + +\item \texttt{yv}---is the y-coordinate for the vertex. + +\item \texttt{zv}---is the z-coordinate for the vertex. + +\end{description} +\item \textbf{Block: CELL2D} + +\begin{description} +\item \texttt{icell2d}---is the cell2d number. Records in the cell2d block must be listed in consecutive order from the first to the last. + +\item \texttt{fdc}---is the fractional distance to the cell center. FDC is relative to the first vertex in the ICVERT array. In most cases FDC should be 0.5, which would place the center of the line segment that defines the cell. If the value of FDC is 1, the cell center would located at the last vertex. FDC values of 0 and 1 can be used to place the node at either end of the cell which can be useful for cells with boundary conditions. + +\item \texttt{ncvert}---is the number of vertices required to define the cell. There may be a different number of vertices for each cell. + +\item \texttt{icvert}---is an array of integer values containing vertex numbers (in the VERTICES block) used to define the cell. Vertices must be listed in the order that defines the line representing the cell. Cells that are connected must share vertices. The bottom elevation of the cell is calculated using the ZV of the first and last vertex point and FDC. + +\end{description} + diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-dimensions.dat b/doc/mf6io/mf6ivar/tex/snf-disl-dimensions.dat new file mode 100644 index 00000000000..f823db64bc0 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-dimensions.dat @@ -0,0 +1,4 @@ +BEGIN DIMENSIONS + NODES + [NVERT ] +END DIMENSIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-griddata.dat b/doc/mf6io/mf6ivar/tex/snf-disl-griddata.dat new file mode 100644 index 00000000000..d25bc8ec829 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-griddata.dat @@ -0,0 +1,8 @@ +BEGIN GRIDDATA + SEGMENT_LENGTH + -- READARRAY + TOSEGMENT + -- READARRAY + [IDOMAIN + -- READARRAY] +END GRIDDATA diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-options.dat b/doc/mf6io/mf6ivar/tex/snf-disl-options.dat new file mode 100644 index 00000000000..7293142029f --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-options.dat @@ -0,0 +1,9 @@ +BEGIN OPTIONS + [LENGTH_UNITS ] + [LENGTH_CONVERT ] + [TIME_CONVERT ] + [NOGRB] + [XORIGIN ] + [YORIGIN ] + [ANGROT ] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-disl-vertices.dat b/doc/mf6io/mf6ivar/tex/snf-disl-vertices.dat new file mode 100644 index 00000000000..76c179add38 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-disl-vertices.dat @@ -0,0 +1,5 @@ +BEGIN VERTICES + + + ... +END VERTICES diff --git a/doc/mf6io/mf6ivar/tex/snf-flw-desc.tex b/doc/mf6io/mf6ivar/tex/snf-flw-desc.tex new file mode 100644 index 00000000000..8010a176bfd --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-flw-desc.tex @@ -0,0 +1,49 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\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{auxmultname}---name of auxiliary variable to be used as multiplier of flow rate. + +\item \texttt{BOUNDNAMES}---keyword to indicate that boundary names may be provided with the list of inflow cells. + +\item \texttt{PRINT\_INPUT}---keyword to indicate that the list of inflow information will be written to the listing file immediately after it is read. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that the list of inflow 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 inflow flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. + +\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 inflow 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 inflow package. + +\end{description} +\item \textbf{Block: DIMENSIONS} + +\begin{description} +\item \texttt{maxbound}---integer value specifying the maximum number of inflow cells that will be specified for use during any stress period. + +\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{cellid}---is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. + +\item \textcolor{blue}{\texttt{q}---is the volumetric inflow rate. A positive value indicates inflow to the stream. Negative values are not allows. 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{aux}---represents the values of the auxiliary variables for each inflow. The values of auxiliary variables must be present for each inflow. 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 inflow 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} + diff --git a/doc/mf6io/mf6ivar/tex/snf-flw-dimensions.dat b/doc/mf6io/mf6ivar/tex/snf-flw-dimensions.dat new file mode 100644 index 00000000000..7b4c7bf6ec7 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-flw-dimensions.dat @@ -0,0 +1,3 @@ +BEGIN DIMENSIONS + MAXBOUND +END DIMENSIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-flw-options.dat b/doc/mf6io/mf6ivar/tex/snf-flw-options.dat new file mode 100644 index 00000000000..0985bd51e40 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-flw-options.dat @@ -0,0 +1,10 @@ +BEGIN OPTIONS + [AUXILIARY ] + [AUXMULTNAME ] + [BOUNDNAMES] + [PRINT_INPUT] + [PRINT_FLOWS] + [SAVE_FLOWS] + [TS6 FILEIN ] + [OBS6 FILEIN ] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-flw-period.dat b/doc/mf6io/mf6ivar/tex/snf-flw-period.dat new file mode 100644 index 00000000000..8ba81439e8c --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-flw-period.dat @@ -0,0 +1,5 @@ +BEGIN PERIOD + <@q@> [<@aux(naux)@>] [] + <@q@> [<@aux(naux)@>] [] + ... +END PERIOD diff --git a/doc/mf6io/mf6ivar/tex/snf-mmr-desc.tex b/doc/mf6io/mf6ivar/tex/snf-mmr-desc.tex new file mode 100644 index 00000000000..e0cb3847ed2 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-mmr-desc.tex @@ -0,0 +1,25 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\item \texttt{SAVE\_FLOWS}---keyword to indicate that budget flow terms will be written to the file specified with ``BUDGET SAVE FILE'' in Output Control. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that calculated flows between cells 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. This option can produce extremely large list files because all cell-by-cell flows are printed. It should only be used with the MMR Package for models that have a small number of cells. + +\end{description} +\item \textbf{Block: GRIDDATA} + +\begin{description} +\item \texttt{iseg\_order}---segment calculation order + +\item \texttt{mann\_n}---manning roughness coefficient + +\item \texttt{seg\_depth}---Depth of bankfull water in segment + +\item \texttt{seg\_slope}---surface slope of each segment + +\item \texttt{x\_coef}---The amount of attenuation of the flow wave, called the Muskingum routing weighting factor; enter 0.0 for reservoirs, diversions, and segment(s) flowing out of the basin + +\end{description} + diff --git a/doc/mf6io/mf6ivar/tex/snf-mmr-griddata.dat b/doc/mf6io/mf6ivar/tex/snf-mmr-griddata.dat new file mode 100644 index 00000000000..db950a2c124 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-mmr-griddata.dat @@ -0,0 +1,12 @@ +BEGIN GRIDDATA + ISEG_ORDER + -- READARRAY + MANN_N + -- READARRAY + SEG_DEPTH + -- READARRAY + SEG_SLOPE + -- READARRAY + X_COEF + -- READARRAY +END GRIDDATA diff --git a/doc/mf6io/mf6ivar/tex/snf-mmr-options.dat b/doc/mf6io/mf6ivar/tex/snf-mmr-options.dat new file mode 100644 index 00000000000..1beb8894206 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-mmr-options.dat @@ -0,0 +1,4 @@ +BEGIN OPTIONS + [SAVE_FLOWS] + [PRINT_FLOWS] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-nam-desc.tex b/doc/mf6io/mf6ivar/tex/snf-nam-desc.tex new file mode 100644 index 00000000000..e629867161a --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-nam-desc.tex @@ -0,0 +1,25 @@ +% DO NOT MODIFY THIS FILE DIRECTLY. IT IS CREATED BY mf6ivar.py + +\item \textbf{Block: OPTIONS} + +\begin{description} +\item \texttt{list}---is name of the listing file to create for this SNF model. If not specified, then the name of the list file will be the basename of the SNF model name file and the '.lst' extension. For example, if the SNF name file is called ``my.model.nam'' then the list file will be called ``my.model.lst''. + +\item \texttt{PRINT\_INPUT}---keyword to indicate that the list of all model stress package information will be written to the listing file immediately after it is read. + +\item \texttt{PRINT\_FLOWS}---keyword to indicate that the list of all model package flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. + +\item \texttt{SAVE\_FLOWS}---keyword to indicate that all model package flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. + +\end{description} +\item \textbf{Block: PACKAGES} + +\begin{description} +\item \texttt{ftype}---is the file type, which must be one of the following character values shown in table~\ref{table:ftype}. Ftype may be entered in any combination of uppercase and lowercase. + +\item \texttt{fname}---is the name of the file containing the package input. The path to the file should be included if the file is not located in the folder where the program was run. + +\item \texttt{pname}---is the user-defined name for the package. PNAME is restricted to 16 characters. No spaces are allowed in PNAME. PNAME character values are read and stored by the program for stress packages only. These names may be useful for labeling purposes when multiple stress packages of the same type are located within a single SNF Model. If PNAME is specified for a stress package, then PNAME will be used in the flow budget table in the listing file; it will also be used for the text entry in the cell-by-cell budget file. PNAME is case insensitive and is stored in all upper case letters. + +\end{description} + diff --git a/doc/mf6io/mf6ivar/tex/snf-nam-options.dat b/doc/mf6io/mf6ivar/tex/snf-nam-options.dat new file mode 100644 index 00000000000..a65ebd5e24d --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-nam-options.dat @@ -0,0 +1,6 @@ +BEGIN OPTIONS + [LIST ] + [PRINT_INPUT] + [PRINT_FLOWS] + [SAVE_FLOWS] +END OPTIONS diff --git a/doc/mf6io/mf6ivar/tex/snf-nam-packages.dat b/doc/mf6io/mf6ivar/tex/snf-nam-packages.dat new file mode 100644 index 00000000000..ee5dc814ee7 --- /dev/null +++ b/doc/mf6io/mf6ivar/tex/snf-nam-packages.dat @@ -0,0 +1,5 @@ +BEGIN PACKAGES + [] + [] + ... +END PACKAGES diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 6d7f29668ac..dff13dc6440 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -131,6 +131,7 @@ + @@ -147,6 +148,14 @@ + + + + + + + + @@ -158,6 +167,7 @@ + @@ -186,6 +196,9 @@ + + + @@ -193,9 +206,6 @@ - - - diff --git a/src/Model/BaseModel.f90 b/src/Model/BaseModel.f90 index 29a71992fd2..ba10ac6d8c1 100644 --- a/src/Model/BaseModel.f90 +++ b/src/Model/BaseModel.f90 @@ -27,7 +27,11 @@ module BaseModelModule procedure :: model_df procedure :: model_ar procedure :: model_rp + procedure :: model_ad procedure :: model_calculate_delt + procedure :: model_solve + procedure :: model_cq + procedure :: model_bd procedure :: model_ot procedure :: model_fp procedure :: model_da @@ -79,6 +83,20 @@ subroutine model_rp(this) return end subroutine model_rp + subroutine model_ad(this) +! ****************************************************************************** +! model_ad -- Advance +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(BaseModelType) :: this +! ------------------------------------------------------------------------------ + ! + ! -- return + return + end subroutine model_ad + subroutine model_calculate_delt(this) ! ****************************************************************************** ! model_calculate_delt -- Calculate time step length @@ -93,6 +111,52 @@ subroutine model_calculate_delt(this) return end subroutine model_calculate_delt + subroutine model_solve(this) +! ****************************************************************************** +! model_solve -- Solve the model (for explicit models) +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(BaseModelType) :: this +! ------------------------------------------------------------------------------ + ! + ! -- return + return + end subroutine model_solve + + subroutine model_cq(this, icnvg, isuppress_output) +! ****************************************************************************** +! model_cq -- Calculate flow +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(BaseModelType) :: this + integer(I4B), intent(in) :: icnvg + integer(I4B), intent(in) :: isuppress_output +! ------------------------------------------------------------------------------ + ! + ! -- return + return + end subroutine model_cq + + subroutine model_bd(this, icnvg, isuppress_output) +! ****************************************************************************** +! model_bd -- Model budget +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(BaseModelType) :: this + integer(I4B), intent(in) :: icnvg + integer(I4B), intent(in) :: isuppress_output +! ------------------------------------------------------------------------------ + ! + ! -- return + return + end subroutine model_bd + subroutine model_ot(this) ! ****************************************************************************** ! model_ot -- output results diff --git a/src/Model/ModelUtilities/Connections.f90 b/src/Model/ModelUtilities/Connections.f90 index de138daabcb..d4219ca82be 100644 --- a/src/Model/ModelUtilities/Connections.f90 +++ b/src/Model/ModelUtilities/Connections.f90 @@ -45,6 +45,7 @@ module ConnectionsModule procedure :: disconnections procedure :: disvconnections procedure :: disuconnections + procedure :: dislconnections procedure :: iajausr procedure :: getjaindex procedure :: set_mask @@ -997,6 +998,98 @@ subroutine disuconnections(this, name_model, nodes, nodesuser, nrsize, & return end subroutine disuconnections + !> @brief procedure to fill the connections object for a disl package + !! + !! todo: No handling yet of cl1, cl2, hwva, etc. for disl as they are not needed + !! and only unreduced disl grids are allowed at the moment + !! + !< + subroutine dislconnections(this, name_model, tosegment) + ! -- modules + use MemoryManagerModule, only: mem_deallocate, mem_setptr + use SparseModule, only: sparsematrix + ! -- dummy + class(ConnectionsType) :: this + character(len=*), intent(in) :: name_model + integer(I4B), dimension(:), intent(in) :: tosegment + ! -- local + type(sparsematrix) :: sparse + integer(I4B) :: ierror + ! + ! -- Allocate scalars + call this%allocate_scalars(name_model) + ! + ! -- Set scalars + this%nodes = size(tosegment) + this%ianglex = 0 + ! + ! -- create sparse matrix object using tosegment + ! and fill ia and ja + call sparse_from_tosegment(tosegment, sparse) + this%nja = sparse%nnz + this%njas = (this%nja - this%nodes) / 2 + ! + ! -- Allocate index arrays of size nja and symmetric arrays + call this%allocate_arrays() + ! + ! -- Fill the IA and JA arrays from sparse, then destroy sparse + call sparse%sort() + call sparse%filliaja(this%ia, this%ja, ierror) + call sparse%destroy() + ! + ! -- fill the isym and jas arrays + call fillisym(this%nodes, this%nja, this%ia, this%ja, this%isym) + call filljas(this%nodes, this%nja, this%ia, this%ja, this%isym, this%jas) + ! + ! -- If reduced system, then need to build iausr and jausr, otherwise point + ! them to ia and ja. + ! TODO: not handled yet for reduced system + !this%iausr => this%ia + !this%jausr => this%ja + !call this%iajausr(nrsize, nodesuser, nodereduced, nodeuser) + ! -- iausr and jausr will be pointers + call mem_deallocate(this%iausr) + call mem_deallocate(this%jausr) + call mem_setptr(this%iausr, 'IA', this%memoryPath) + call mem_setptr(this%jausr, 'JA', this%memoryPath) + + return + end subroutine dislconnections + + !> @brief Using tosegment, fill the sparse object + !< + subroutine sparse_from_tosegment(tosegment, spm) + ! -- modules + use SparseModule, only: sparsematrix + ! -- dummy + integer(I4B), dimension(:), contiguous :: tosegment + type(sparsematrix), intent(inout) :: spm + ! -- local + integer(I4B) :: nodes + integer(I4B) :: n + integer(I4B) :: j + + nodes = size(tosegment) + call spm%init(nodes, nodes, 3) + ! + ! -- insert diagonal so it stays in front + do n = 1, size(tosegment) + call spm%addconnection(n, n, 1) + end do + ! + ! -- add tosegment connections for non-zero + ! segments (0 indicates no downstream segment) + do n = 1, size(tosegment) + j = tosegment(n) + if (j > 0) then + call spm%addconnection(n, j, 1) + call spm%addconnection(j, n, 1) + end if + end do + + return + end subroutine sparse_from_tosegment + subroutine iajausr(this, nrsize, nodesuser, nodereduced, nodeuser) ! ****************************************************************************** ! iajausr -- Fill iausr and jausr if reduced grid, otherwise point them diff --git a/src/Model/ModelUtilities/DislGeom.f90 b/src/Model/ModelUtilities/DislGeom.f90 new file mode 100644 index 00000000000..5f6ca167580 --- /dev/null +++ b/src/Model/ModelUtilities/DislGeom.f90 @@ -0,0 +1,416 @@ +module DislGeom + + use KindModule, only: DP, I4B + use SimModule, only: store_error, store_error_unit, ustop + implicit none + private + public :: DislGeomType + public :: calcdist, partialdist, line_unit_vector + + type DislGeomType + integer(I4B) :: nodesuser + logical :: reduced + integer(I4B) :: nodes ! number of reduced nodes; + real(DP), pointer, dimension(:) :: cellfdc => null() + integer(I4B), pointer, dimension(:) :: iavert => null() + integer(I4B), pointer, dimension(:) :: javert => null() + integer(I4B), pointer, dimension(:) :: iavertcells => null() + integer(I4B), pointer, dimension(:) :: javertcells => null() + real(DP), pointer, dimension(:, :) :: vertices => null() + real(DP), pointer, dimension(:, :) :: cellcenters => null() + integer(I4B), pointer, dimension(:, :) :: centerverts => null() + integer(I4B), pointer, dimension(:) :: nodereduced => null() + integer(I4B), pointer, dimension(:) :: nodeuser => null() + contains + procedure :: init + procedure :: cprops + procedure :: connection_vector + procedure :: shared_vertex + procedure :: disttocenter + procedure :: containscenter + procedure :: vertccdist + end type DislGeomType + + contains + + subroutine init(this, nodesuser, nodes, cellfdc, iavert, javert, & + iavertcells, javertcells, vertices, cellcenters, & + centerverts, nodereduced, nodeuser) + class(DislGeomType) :: this + integer(I4B), intent(in) :: nodesuser + integer(I4B), intent(in) :: nodes + real(DP), dimension(nodes), target :: cellfdc + integer(I4B), dimension(:), target :: iavert + integer(I4B), dimension(:), target :: javert + integer(I4B), dimension(:), target :: iavertcells + integer(I4B), dimension(:), target :: javertcells + real(DP), dimension(:, :), target :: vertices + real(DP), dimension(:, :), target :: cellcenters + integer(I4B), dimension(:, :), target :: centerverts + integer(I4B), dimension(:), target :: nodereduced + integer(I4B), dimension(:), target :: nodeuser + ! -- local + this%nodes = nodes + this%nodesuser = nodesuser + this%cellfdc => cellfdc + this%iavert => iavert + this%javert => javert + this%iavertcells => iavertcells + this%javertcells => javertcells + + this%vertices => vertices + this%cellcenters => cellcenters + this%centerverts => centerverts + this%nodereduced => nodereduced + this%nodeuser => nodeuser + if(nodes < nodesuser) then + this%reduced = .true. + else + this%reduced = .false. + endif + end subroutine init + + subroutine cprops(this, cell1, cell2, hwva, cl1, cl2) + ! -- module + use ConstantsModule, only: DZERO, DHALF, DONE + class(DislGeomType) :: this + integer(I4B), intent(in) :: cell1 + integer(I4B), intent(in) :: cell2 + real(DP), intent(out) :: hwva + real(DP), intent(out) :: cl1 + real(DP), intent(out) :: cl2 + ! -- local + integer(I4B) :: shared_vert + character(len=300) :: ermsg + character(len=*), parameter :: fmtvert = & + "('ERROR. CELLS ', I0, ' AND ', i0, ' DO NOT SHARE A VERTEX ')" + + ! find shared vertex + call this%shared_vertex(cell1, cell2, shared_vert) + if (shared_vert == 0) then + write(ermsg, fmtvert) cell1, cell2 + call store_error(ermsg) + call ustop() + end if + + ! cl1 - distance from center of pipe 1 to vertex + cl1 = this%disttocenter(cell1, shared_vert) + + ! cl2 - distance from center of pipe 2 to connection + cl2 = this%disttocenter(cell2, shared_vert) + + ! hwva - smallest cross sectional area? + + return + end subroutine cprops + + subroutine shared_vertex(this, cell1, cell2, shared_vert) + ! return the shared vertex of two cells + ! -- module + class(DislGeomType) :: this + integer(I4B), intent(in) :: cell1 + integer(I4B), intent(in) :: cell2 + integer(I4B), intent(out) :: shared_vert + ! -- local + integer(I4B) :: n, v, test1, test3, test4, test5 + ! + ! find shared vertex + shared_vert = 0 + ! loop through all vertices in cell1 + outer: do v = this%iavert(cell1), this%iavert(cell1+1) - 1 + ! loop through all cells that share vertex v + test1 = this%javert(v) + test3 = this%iavertcells(this%javert(v)) + test4 = this%iavertcells(this%javert(v)+1) + do n = this%iavertcells(this%javert(v)), this%iavertcells(this%javert(v)+1) - 1 + ! if cell2 has shared vertex v + test5 = this%javertcells(n) + if (cell2 == this%javertcells(n)) then + ! save shared vertex and exit + shared_vert = this%javert(v) + exit outer + end if + end do + end do outer + return + end subroutine shared_vertex + + subroutine connection_vector(this, cell1, cell2, nozee, xcomp, & + ycomp, zcomp, conlen) + ! return the x y and z components of a unit vector that points + ! from the center of this to the center of cell2, and the + ! straight-line connection length + ! -- module + use ConstantsModule, only: DZERO, DHALF, DONE + ! -- dummy + class(DislGeomType) :: this + integer(I4B), intent(in) :: cell1 + integer(I4B), intent(in) :: cell2 + logical, intent(in) :: nozee + real(DP), intent(out) :: xcomp + real(DP), intent(out) :: ycomp + real(DP), intent(out) :: zcomp + real(DP), intent(out) :: conlen + ! -- local + real(DP) :: x1, y1, z1, x2, y2, z2 + ! + x1 = this%cellcenters(1, cell1) + y1 = this%cellcenters(2, cell1) + x2 = this%cellcenters(1, cell2) + y2 = this%cellcenters(2, cell2) + if (nozee) then + z1 = DZERO + z2 = DZERO + else + z1 = this%cellcenters(3, cell1) + z2 = this%cellcenters(3, cell2) + end if + call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, & + conlen) + return + end subroutine connection_vector + + function disttocenter(this, nodenum, vertnum) result(dist) +! ****************************************************************************** +! disttocenter -- Return distance from cell vertex to cell center +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(DislGeomType) :: this + integer(I4B), intent(in) :: nodenum + integer(I4B), intent(in) :: vertnum + real(DP) :: dist + integer(I4B) :: j + logical :: cellcenter, vertex + character(len=300) :: ermsg + character(len=*), parameter :: novert = & + "('ERROR. VERTEX ', I0, ' NOT FOUND IN NODE ', i0, '.')" + character(len=*), parameter :: cellcent = & + "('ERROR. CELL CENTER NOT FOUND FOR NODE ', I0, '.')" + +! ------------------------------------------------------------------------------ + ! initialize + j = this%iavert(nodenum) + if (this%javert(j) == vertnum) then + vertex = .true. + if (this%containscenter(nodenum, this%javert(j), & + this%javert(j+1)) .eqv. .true.) then + cellcenter = .true. + dist = this%vertccdist(nodenum, this%javert(j)) + return + else + cellcenter = .false. + dist = calcdist(this%vertices, this%javert(j), this%javert(j+1)) + end if + else + vertex = .false. + if (this%containscenter(nodenum, this%javert(j), & + this%javert(j+1)) .eqv. .true.) then + cellcenter = .true. + dist = this%vertccdist(nodenum, this%javert(j+1)) + else + cellcenter = .false. + dist = 0.0 + end if + end if + + ! loop through cell's vertices + loop: do j = this%iavert(nodenum)+1, this%iavert(nodenum+1) - 1 + if (this%javert(j) == vertnum) then + vertex = .true. + if (cellcenter .eqv. .true.) then + exit loop + else + if (this%containscenter(nodenum, this%javert(j), & + this%javert(j+1)) .eqv. .true.) then + cellcenter = .true. + ! add the distance from this vertex to the cell center + dist = dist + this%vertccdist(nodenum, this%javert(j)) + exit loop + else + ! add the distance between this vertex and the next + dist = dist + calcdist(this%vertices, this%javert(j), & + this%javert(j+1)) + end if + end if + else if (this%containscenter(nodenum, this%javert(j), & + this%javert(j+1)) .eqv. .true.) then + cellcenter = .true. + if (vertex .eqv. .true.) then + ! add the remaining distance from the current vertex to the cell center + dist = dist + this%vertccdist(nodenum, this%javert(j)) + exit loop + else + ! add the distance from the cell center to the next vertex + dist = dist + this%vertccdist(nodenum, this%javert(j+1)) + end if + else if ((vertex .eqv. .true.) .or. (cellcenter .eqv. .true.)) then + ! add the distance between this vertex and the next + dist = dist + calcdist(this%vertices, this%javert(j), & + this%javert(j+1)) + end if + end do loop + + if (vertex .eqv. .false.) then + write(ermsg, novert) vertnum, nodenum + call store_error(ermsg) + call ustop() + else if (cellcenter .eqv. .false.) then + write(ermsg, cellcent) nodenum + call store_error(ermsg) + call ustop() + end if + + return + end function disttocenter + + function containscenter(this, nodenum, vert1, vert2) result(l) +! ****************************************************************************** +! containscenter -- check to see if cell center is contained at or +! between vert1 and vert2 + +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + class(DislGeomType) :: this + integer(I4B), intent(in) :: nodenum + integer(I4B), intent(in) :: vert1 + integer(I4B), intent(in) :: vert2 + logical :: l +! ------------------------------------------------------------------------------ + if (this%centerverts(2, nodenum) == 0) then + ! cell center is at a vertex + if ((this%centerverts(1, nodenum) == vert1) .or. & + (this%centerverts(1, nodenum) == vert2)) then + l = .true. + else + l = .false. + end if + else + ! cell center is between two vertices + if (((this%centerverts(1, nodenum) == vert1) .or. & + (this%centerverts(1, nodenum) == vert2)) .and. & + ((this%centerverts(2, nodenum) == vert1) .or. & + (this%centerverts(2, nodenum) == vert2))) then + l = .true. + else + l = .false. + end if + end if + + return + end function containscenter + + function vertccdist(this, node, vertex) result(dist) +! ****************************************************************************** +! vertccdist -- Return the distance between a cell vertex and that cell's center. +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + ! -- dummy + class(DislGeomType), intent(in) :: this + integer(I4B), intent(in) :: node + integer(I4B), intent(in) :: vertex + real(DP) :: dist + real(DP) :: xdist, ydist, zdist + ! -- local +! ------------------------------------------------------------------------------ + ! + ! -- calc distance + xdist = abs(this%vertices(1, vertex) - this%cellcenters(1, node)) + ydist = abs(this%vertices(2, vertex) - this%cellcenters(2, node)) + zdist = abs(this%vertices(3, vertex) - this%cellcenters(3, node)) + dist = sqrt(xdist * xdist + ydist * ydist + zdist * zdist) + + ! -- return + return + end function vertccdist + + function calcdist(vertices, vert1, vert2) result(dist) +! ****************************************************************************** +! calcdist -- Return the distance between two vertices. +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + ! -- dummy + real(DP), dimension(:, :), intent(in) :: vertices + integer(I4B), intent(in) :: vert1 + integer(I4B), intent(in) :: vert2 + real(DP) :: dist, xdist, ydist, zdist + ! -- local +! ------------------------------------------------------------------------------ + ! + ! -- calc distance + xdist = abs(vertices(1, vert1) - vertices(1, vert2)) + ydist = abs(vertices(2, vert1) - vertices(2, vert2)) + zdist = abs(vertices(3, vert1) - vertices(3, vert2)) + dist = sqrt(xdist * xdist + ydist * ydist + zdist * zdist) + + ! -- return + return + end function calcdist + + function partialdist(coord1, coord2, percent) result(dist) +! ****************************************************************************** +! partialdist -- Return the distance between two vertices. +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + use ConstantsModule, only: DONE + ! -- dummy + real(DP), intent(in) :: coord1 + real(DP), intent(in) :: coord2 + real(DP) :: percent + real(DP) :: dist + ! -- local +! ------------------------------------------------------------------------------ + ! + ! -- calc distance + if (coord1 > coord2) then + dist = coord2 + (coord1 - coord2) * (DONE - percent) + else + dist = coord1 + (coord2 - coord1) * percent + end if + + ! -- return + return + end function partialdist + + subroutine line_unit_vector(x0, y0, z0, x1, y1, z1, & + xcomp, ycomp, zcomp, vmag) +! ****************************************************************************** +! line_unit_vector -- Calculate the vector components (xcomp, ycomp, and zcomp) +! for a line defined by two points, (x0, y0, z0), (x1, y1, z1). Also return +! the magnitude of the original vector, vmag. +! ****************************************************************************** +! +! SPECIFICATIONS: +! ------------------------------------------------------------------------------ + real(DP), intent(in) :: x0 + real(DP), intent(in) :: y0 + real(DP), intent(in) :: z0 + real(DP), intent(in) :: x1 + real(DP), intent(in) :: y1 + real(DP), intent(in) :: z1 + real(DP), intent(out) :: xcomp + real(DP), intent(out) :: ycomp + real(DP), intent(out) :: zcomp + real(DP) :: dx, dy, dz, vmag +! ------------------------------------------------------------------------------ + dx = x1 - x0 + dy = y1 - y0 + dz = z1 - z0 + vmag = sqrt(dx ** 2 + dy ** 2 + dz ** 2) + xcomp = dx / vmag + ycomp = dy / vmag + zcomp = dz / vmag + return + end subroutine line_unit_vector + + +end module DislGeom \ No newline at end of file diff --git a/src/Model/NumericalModel.f90 b/src/Model/NumericalModel.f90 index 96b10b221aa..cec5378235b 100644 --- a/src/Model/NumericalModel.f90 +++ b/src/Model/NumericalModel.f90 @@ -468,6 +468,8 @@ subroutine load_input_context(this, filtyp, modelname, pkgname, inunit, iout, & call input_load('DISU6', 'GWF', 'DISU', modelname, pkgname, inunit, iout) case ('DISV6') call input_load('DISV6', 'GWF', 'DISV', modelname, pkgname, inunit, iout) + case ('DISL6') + call input_load('DISL6', 'SNF', 'DISL', modelname, pkgname, inunit, iout) case default end select ! diff --git a/src/Model/StreamNetworkFlow/snf1.f90 b/src/Model/StreamNetworkFlow/snf1.f90 new file mode 100644 index 00000000000..4cbec10f5af --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1.f90 @@ -0,0 +1,900 @@ + +!> @brief Stream Network Flow (SNF) Module +!! +!! This module contains the SNF Model +!! +!! Status and remaining tasks +!! ONGOING -- Implement SNF infrastructure +!! DONE -- Implement Explicit Model Solution (EMS6) to handle explicit models +!! DONE -- Implement DISL Package +!! ONGOING -- Implement MMR Package +!! DONE -- Use MMR package to solve a flow problem in combination with DISL and FLW +!! DONE -- Implement FLW Package to handle lateral and point inflows +!! DONE -- Transfer results into the flowja vector +!! DONE -- Implement strategy for storing outflow terms and getting them into budget +!! DONE -- Implement SNF and FLW advance routines to handle transient problems +!! DONE -- Implement storage terms and getting them into budget +!! DONE -- Observations +!! DONE -- Initial conditions? +!! DONE -- Rework the Iterative Model Solution (IMS6) to handle both implicit and explicit models +! Look into mass conservative MC method (https://hess.copernicus.org/articles/11/1645/2007/hess-11-1645-2007.pdf) +!! Deal with the timestep and subtiming issues +!! Implement output control +!! Flopy support for DISL and DISL binary grid file +!! mf6io guide +!! Mover support? +!! SNF-SNF Exchange +!! SNF-SNF Exchange in parallel +!! Use dag_module to calculate iseg_order (if iseg_order not specified by user) +!! +!< +module SnfModule + + use KindModule, only: DP, I4B + use ConstantsModule, only: DZERO, LENFTYPE, DNODATA + use InputOutputModule, only: ParseLine, upcase + use SimVariablesModule, only: errmsg + use MemoryManagerModule, only: mem_allocate + use VersionModule, only: write_listfile_header + use BaseModelModule, only: BaseModelType + use NumericalModelModule, only: NumericalModelType + use BndModule, only: BndType, AddBndToList, GetBndFromList + use SnfMmrModule, only: SnfMmrType + use SnfOcModule, only: SnfOcType + use BudgetModule, only: BudgetType + + implicit none + + private + public :: snf_cr + + type, extends(NumericalModelType) :: SnfModelType + type(SnfMmrType), pointer :: mmr => null() !< muskingum routing package + type(SnfOcType), pointer :: oc => null() ! output control package + type(BudgetType), pointer :: budget => null() ! budget object + integer(I4B), pointer :: inmmr => null() ! unit number MMR + integer(I4B), pointer :: inoc => null() ! unit number OC + contains + procedure :: allocate_scalars + procedure :: allocate_arrays + procedure :: model_df => snf_df + procedure :: model_ar => snf_ar + procedure :: model_rp => snf_rp + procedure :: model_ad => snf_ad + procedure :: model_solve => snf_solve !< routine for solving this model for the current time step + procedure :: model_cq => snf_cq + procedure :: model_bd => snf_bd + procedure :: model_ot => snf_ot + procedure :: model_da => snf_da + procedure :: snf_ot_obs + procedure :: snf_ot_flow + procedure :: snf_ot_dv + procedure :: snf_ot_bdsummary + procedure :: package_create + procedure :: ftype_check + procedure :: load_input_context => snf_load_input_context + end type SnfModelType + + ! -- Module variables constant for simulation + integer(I4B), parameter :: NIUNIT = 100 + character(len=LENFTYPE), dimension(NIUNIT) :: cunit + data cunit/'DISL6', 'MMR6 ', 'OC6 ', 'FLW6 ', ' ', & ! 5 + &95*' '/ + + contains + + !> @brief Create a new stream network flow model object + !! + !! (1) creates model object and add to modellist + !! (2) assign values + !! + !< + subroutine snf_cr(filename, id, modelname) + ! -- modules + use ListsModule, only: basemodellist + use MemoryHelperModule, only: create_mem_path + use BaseModelModule, only: AddBaseModelToList + use SimModule, only: store_error, count_errors + use GenericUtilitiesModule, only: write_centered + use ConstantsModule, only: LINELENGTH, LENPACKAGENAME + use MemoryManagerModule, only: mem_allocate + use SnfDislModule, only: disl_cr + use SnfMmrModule, only: mmr_cr + use SnfOcModule, only: oc_cr + use BudgetModule, only: budget_cr + use NameFileModule, only: NameFileType + ! -- dummy + character(len=*), intent(in) :: filename + integer(I4B), intent(in) :: id + character(len=*), intent(in) :: modelname + ! -- local + integer(I4B) :: indisl + integer(I4B) :: ipakid, i, j, iu, ipaknum + character(len=LENPACKAGENAME) :: pakname + type(NameFileType) :: namefile_obj + type(SnfModelType), pointer :: this + class(BaseModelType), pointer :: model + integer(I4B) :: nwords + character(len=LINELENGTH), allocatable, dimension(:) :: words + ! -- format + ! + ! -- Allocate a new SNF Model (this) and add it to basemodellist + allocate (this) + ! + ! -- Set memory path before allocation in memory manager can be done + this%memoryPath = create_mem_path(modelname) + ! + call this%allocate_scalars(modelname) + model => this + call AddBaseModelToList(basemodellist, model) + ! + ! -- Assign values + this%filename = filename + this%name = modelname + this%macronym = 'SNF' + this%id = id + ! + ! -- Open namefile and set iout + call namefile_obj%init(this%filename, 0) + call namefile_obj%add_cunit(niunit, cunit) + call namefile_obj%openlistfile(this%iout) + ! + ! -- Write header to model list file + call write_listfile_header(this%iout, 'STREAM NETWORK FLOW MODEL (SNF)') + ! + ! -- Open files + call namefile_obj%openfiles(this%iout) + ! + ! -- SNF options + if (size(namefile_obj%opts) > 0) then + write (this%iout, '(1x,a)') 'NAMEFILE OPTIONS:' + end if + ! + ! -- Parse options in the SNF name file + do i = 1, size(namefile_obj%opts) + call ParseLine(namefile_obj%opts(i), nwords, words) + call upcase(words(1)) + select case (words(1)) + case ('PRINT_INPUT') + this%iprpak = 1 + write (this%iout, '(4x,a)') 'STRESS PACKAGE INPUT WILL BE PRINTED '// & + 'FOR ALL MODEL STRESS PACKAGES' + case ('PRINT_FLOWS') + this%iprflow = 1 + write (this%iout, '(4x,a)') 'PACKAGE FLOWS WILL BE PRINTED '// & + 'FOR ALL MODEL PACKAGES' + case ('SAVE_FLOWS') + this%ipakcb = -1 + write (this%iout, '(4x,a)') & + 'FLOWS WILL BE SAVED TO BUDGET FILE SPECIFIED IN OUTPUT CONTROL' + case default + write (errmsg, '(4x,a,a,a,a)') & + 'Unknown SNF namefile (', & + trim(adjustl(this%filename)), ') option: ', & + trim(adjustl(namefile_obj%opts(i))) + call store_error(errmsg, terminate=.TRUE.) + end select + end do + ! + ! -- Assign unit numbers to attached modules, and remove + ! -- from unitnumber (by specifying 1 for iremove) + ! + indisl = 0 + call namefile_obj%get_unitnumber('DISL6', indisl, 1) + call namefile_obj%get_unitnumber('MMR6', this%inmmr, 1) + call namefile_obj%get_unitnumber('OC6', this%inoc, 1) + ! call namefile_obj%get_unitnumber('OBS6', this%inobs, 1) + ! + ! -- Check to make sure that required ftype's have been specified + call this%ftype_check(namefile_obj, indisl) + ! + ! -- Create discretization object + if (indisl > 0) then + call this%load_input_context('DISL6', this%name, 'DISL', indisl, & + this%iout) + call disl_cr(this%dis, this%name, indisl, this%iout) + end if + ! + ! -- Create utility objects + call budget_cr(this%budget, this%name) + ! + ! -- Muskingum Manning Routing Package + if (this%inmmr > 0) then + call this%load_input_context('MMR6', this%name, 'MMR', this%inmmr, this%iout) + call mmr_cr(this%mmr, this%name, this%inmmr, this%dis, this%iout) + end if + ! + ! -- Create packages that are tied directly to model + ! call ic_cr(this%ic, this%name, this%inic, this%iout, this%dis) + call oc_cr(this%oc, this%name, this%inoc, this%iout) + ! call snf_obs_cr(this%obs, this%inobs) + ! + ! -- Create stress packages + ipakid = 1 + do i = 1, niunit + ipaknum = 1 + do j = 1, namefile_obj%get_nval_for_row(i) + iu = namefile_obj%get_unitnumber_rowcol(i, j) + call namefile_obj%get_pakname(i, j, pakname) + call this%package_create(cunit(i), ipakid, ipaknum, pakname, iu, & + this%iout) + ipaknum = ipaknum + 1 + ipakid = ipakid + 1 + end do + end do + + ! put this here in CR, but it is supposed to be in df. It seems + ! like now that dis can be loaded entirely up front in CR, it may + ! be okay to leave this here now. + ! + ! -- Assign or point model members to dis members + ! this%neq will be incremented if packages add additional unknowns + this%neq = this%dis%nodes + this%nja = this%dis%nja + this%ia => this%dis%con%ia + this%ja => this%dis%con%ja + ! + ! -- Allocate model arrays, now that neq and nja are known + call this%allocate_arrays() + ! + ! -- return + return + end subroutine snf_cr + + !> @brief Allocate memory for scalar members + subroutine allocate_scalars(this, modelname) + ! -- modules + ! -- dummy + class(SnfModelType) :: this + character(len=*), intent(in) :: modelname + ! + ! -- allocate members from parent class + call this%NumericalModelType%allocate_scalars(modelname) + ! + ! -- allocate members that are part of model class + call mem_allocate(this%inmmr, 'INMMR', this%memoryPath) + call mem_allocate(this%inoc, 'INOC', this%memoryPath) + ! + this%inmmr = 0 + this%inoc = 0 + ! + ! -- return + return + end subroutine allocate_scalars + + !> @brief Allocate memory for scalar members + subroutine allocate_arrays(this) + ! -- modules + ! -- dummy + class(SnfModelType) :: this + integer(I4B) :: i + ! + ! -- allocate members from parent class + call this%NumericalModelType%allocate_arrays() + ! + ! -- This is not a numerical solution, so x, rhs, and active + ! are allocated by a numerical solution, so need to do it + ! here. + call mem_allocate(this%x, this%neq, 'X', this%memoryPath) + call mem_allocate(this%rhs, this%neq, 'RHS', this%memoryPath) + call mem_allocate(this%ibound, this%neq, 'IBOUND', this%memoryPath) + ! + ! initialize arrays + do i = 1, this%neq + this%x(i) = DZERO + this%rhs(i) = DZERO + this%ibound(i) = 1 + end do + ! + ! -- return + return + end subroutine allocate_arrays + + !> @brief Define packages of the model + !< + subroutine snf_df(this) + ! -- modules + ! -- dummy + class(SnfModelType) :: this + ! -- local + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! + call this%oc%oc_df() + call this%budget%budget_df(niunit, 'VOLUME', 'L**3') + ! + ! -- Define packages and assign iout for time series managers + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_df(this%neq, this%dis) + end do + ! + ! -- Store information needed for observations + !call this%obs%obs_df(this%iout, this%name, 'SNF', this%dis) + ! + ! -- return + return + end subroutine snf_df + + !> @brief SNF Allocate and Read + !< + subroutine snf_ar(this) + ! -- dummy + class(SnfModelType) :: this + ! -- locals + integer(I4B), dimension(:), allocatable :: itemp + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! -- Allocate and read modules attached to model + ! if (this%inobs > 0) call this%obs%snf_obs_ar(this%ic, this%x, this%flowja) + ! + ! -- Call dis_ar to write binary grid file + allocate(itemp(this%dis%nodes)) + call this%dis%dis_ar(itemp) + call this%mmr%mmr_ar() + deallocate(itemp) + ! + ! -- set up output control + call this%oc%oc_ar(this%x, this%dis, DNODATA) + call this%budget%set_ibudcsv(this%oc%ibudcsv) + ! + ! -- Package input files now open, so allocate and read + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%set_pointers(this%dis%nodes, this%ibound, this%x, & + this%xold, this%flowja) + ! -- Read and allocate package + call packobj%bnd_ar() + end do + ! + ! -- return + return + end subroutine snf_ar + + !> @brief Stream Network Flow Model Read and Prepare + !! + !! (1) calls package read and prepare routines + !! + !< + subroutine snf_rp(this) + ! -- modules + use TdisModule, only: readnewdata + ! -- dummy + class(SnfModelType) :: this + ! -- local + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! + ! -- Check with TDIS on whether or not it is time to RP + if (.not. readnewdata) return + ! + ! -- Read and prepare + if (this%inmmr > 0) call this%mmr%mmr_rp() + if (this%inoc > 0) call this%oc%oc_rp() + !if (this%insto > 0) call this%sto%sto_rp() + !if (this%inmvr > 0) call this%mvr%mvr_rp() + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_rp() + !call packobj%bnd_rp_obs() + end do + ! + ! -- Return + return + end subroutine snf_rp + + !> @brief Stream Network Flow Model Time Step Advance + !! + !! (1) calls package advance subroutines + !! + !< + subroutine snf_ad(this) + ! -- modules + use SimVariablesModule, only: isimcheck, iFailedStepRetry + ! -- dummy + class(SnfModelType) :: this + class(BndType), pointer :: packobj + ! -- local + integer(I4B) :: irestore + integer(I4B) :: ip, n + ! + ! -- Reset state variable + irestore = 0 + if (iFailedStepRetry > 0) irestore = 1 + if (irestore == 0) then + ! + ! -- copy x into xold + do n = 1, this%dis%nodes + this%xold(n) = this%x(n) + end do + else + ! + ! -- copy xold into x if this time step is a redo + do n = 1, this%dis%nodes + this%x(n) = this%xold(n) + end do + end if + ! + ! -- Advance + if (this%inmmr > 0) call this%mmr%mmr_ad(irestore) + !if (this%insto > 0) call this%sto%sto_ad() + !if (this%inmvr > 0) call this%mvr%mvr_ad() + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ad() + if (isimcheck > 0) then + call packobj%bnd_ck() + end if + end do + ! + ! -- Push simulated values to preceding time/subtime step + !call this%obs%obs_ad() + ! + ! -- return + return + end subroutine snf_ad + + !> @brief Make explicit solve for this time step + subroutine snf_solve(this) + ! -- modules + ! -- dummy + class(SnfModelType) :: this + ! -- local + integer(I4B) :: n + integer(I4B) :: i + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! -- Initialize rhs accumulator + do n = 1, this%dis%nodes + this%rhs(n) = DZERO + end do + ! + ! -- Call boundary packages to set up inflows + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_cf() + + ! accumulate individual package rhs into this%rhs + do i = 1, packobj%nbound + n = packobj%nodelist(i) + if (this%ibound(n) > 0) then + this%rhs(n) = this%rhs(n) + packobj%rhs(i) + end if + end do + + end do + + ! -- MMR Solve + call this%mmr%mmr_solve(this%rhs) + + ! -- return + return + end subroutine snf_solve + + !> @brief Calculate flow + !< + subroutine snf_cq(this, icnvg, isuppress_output) + ! -- modules + ! -- dummy + class(SnfModelType) :: this + integer(I4B), intent(in) :: icnvg + integer(I4B), intent(in) :: isuppress_output + ! -- local + integer(I4B) :: i + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! -- Construct the flowja array. Flowja is calculated each time, even if + ! output is suppressed. (flowja is positive into a cell.) The diagonal + ! position of the flowja array will contain the flow residual after + ! these routines are called, so each package is responsible for adding + ! its flow to this diagonal position. + do i = 1, this%nja + this%flowja(i) = DZERO + end do + if (this%inmmr > 0) call this%mmr%mmr_cq(this%flowja) + ! + ! -- Go through packages and call cq routines. cf() routines are called + ! first to regenerate non-linear terms to be consistent with the final + ! head solution. + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_cf(reset_mover=.false.) + call packobj%bnd_cq(this%x, this%flowja) + end do + ! + ! -- Return + return + end subroutine snf_cq + + !> @brief Model Budget + !< + subroutine snf_bd(this, icnvg, isuppress_output) + ! -- modules + use SparseModule, only: csr_diagsum + ! -- dummy + class(SnfModelType) :: this + integer(I4B), intent(in) :: icnvg + integer(I4B), intent(in) :: isuppress_output + ! -- local + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! -- Finalize calculation of flowja by adding face flows to the diagonal. + ! This results in the flow residual being stored in the diagonal + ! position for each cell. + call csr_diagsum(this%dis%con%ia, this%flowja) + ! + ! -- Save the solution convergence flag + this%icnvg = icnvg + ! + ! -- Budget routines (start by resetting). Sole purpose of this section + ! is to add in and outs to model budget. All ins and out for a model + ! should be added here to this%budget. In a subsequent exchange call, + ! exchange flows might also be added. + call this%budget%reset() + if (this%inmmr > 0) call this%mmr%mmr_bd(isuppress_output, this%budget) + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_bd(this%budget) + end do + ! + ! -- Return + return + end subroutine snf_bd + + !> @brief Stream Network Flow Model Output + subroutine snf_ot(this) + ! -- modules + use TdisModule, only: kstp, kper, tdis_ot, endofperiod + ! -- dummy + class(SnfModelType) :: this + ! -- local + integer(I4B) :: idvsave + integer(I4B) :: idvprint + integer(I4B) :: icbcfl + integer(I4B) :: icbcun + integer(I4B) :: ibudfl + integer(I4B) :: ipflag + ! -- formats + character(len=*), parameter :: fmtnocnvg = & + "(1X,/9X,'****FAILED TO MEET SOLVER CONVERGENCE CRITERIA IN TIME STEP ', & + &I0,' OF STRESS PERIOD ',I0,'****')" + ! + ! -- Set write and print flags + idvsave = 0 + idvprint = 0 + icbcfl = 0 + ibudfl = 0 + ! if (this%oc%oc_save('STAGE')) idvsave = 1 + ! if (this%oc%oc_print('STAGE')) idvprint = 1 + if (this%oc%oc_save('BUDGET')) icbcfl = 1 + if (this%oc%oc_print('BUDGET')) ibudfl = 1 + icbcun = this%oc%oc_save_unit('BUDGET') + ! + ! -- Override ibudfl and idvprint flags for nonconvergence + ! and end of period + ibudfl = this%oc%set_print_flag('BUDGET', this%icnvg, endofperiod) + ! idvprint = this%oc%set_print_flag('STAGE', this%icnvg, endofperiod) + ! + ! Calculate and save observations + call this%snf_ot_obs() + ! ! + ! Save and print flows + call this%snf_ot_flow(icbcfl, ibudfl, icbcun) + ! + ! Save and print dependent variables + call this%snf_ot_dv(idvsave, idvprint, ipflag) + ! + ! Print budget summaries + call this%snf_ot_bdsummary(ibudfl, ipflag) + ! + ! -- Timing Output; if any dependendent variables or budgets + ! are printed, then ipflag is set to 1. + if (ipflag == 1) call tdis_ot(this%iout) + ! + ! -- Write non-convergence message + if (this%icnvg == 0) then + write (this%iout, fmtnocnvg) kstp, kper + end if + ! + ! -- Return + return + end subroutine snf_ot + + subroutine snf_ot_obs(this) + class(SnfModelType) :: this + class(BndType), pointer :: packobj + integer(I4B) :: ip + + ! -- Calculate and save mmr observations + if (this%inmmr > 0) then + call this%mmr%mmr_bd_obs() + call this%mmr%obs%obs_ot() + end if + + ! -- Calculate and save package observations + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_bd_obs() + call packobj%bnd_ot_obs() + end do + + end subroutine snf_ot_obs + + subroutine snf_ot_flow(this, icbcfl, ibudfl, icbcun) + class(SnfModelType) :: this + integer(I4B), intent(in) :: icbcfl + integer(I4B), intent(in) :: ibudfl + integer(I4B), intent(in) :: icbcun + class(BndType), pointer :: packobj + integer(I4B) :: ip + + ! -- Save SNF flows + if (this%inmmr > 0) then + call this%mmr%mmr_save_model_flows(this%flowja, icbcfl, icbcun) + end if + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun) + end do + + ! -- Save advanced package flows + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_package_flows(icbcfl=icbcfl, ibudfl=0) + end do + ! if (this%inmvr > 0) then + ! call this%mvr%mvr_ot_saveflow(icbcfl, ibudfl) + ! end if + + ! -- Print SNF flows + if (this%inmmr > 0) call this%mmr%mmr_print_model_flows(ibudfl, this%flowja) + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0) + end do + + ! -- Print advanced package flows + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_package_flows(icbcfl=0, ibudfl=ibudfl) + end do + ! if (this%inmvr > 0) then + ! call this%mvr%mvr_ot_printflow(icbcfl, ibudfl) + ! end if + + end subroutine snf_ot_flow + + subroutine snf_ot_dv(this, idvsave, idvprint, ipflag) + class(SnfModelType) :: this + integer(I4B), intent(in) :: idvsave + integer(I4B), intent(in) :: idvprint + integer(I4B), intent(inout) :: ipflag + class(BndType), pointer :: packobj + integer(I4B) :: ip + ! + ! -- Print advanced package dependent variables + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_dv(idvsave, idvprint) + end do + ! + ! -- save stage and print stage (if implemented) + call this%oc%oc_ot(ipflag) + ! + ! -- Return + return + end subroutine snf_ot_dv + + subroutine snf_ot_bdsummary(this, ibudfl, ipflag) + use TdisModule, only: kstp, kper, totim + class(SnfModelType) :: this + integer(I4B), intent(in) :: ibudfl + integer(I4B), intent(inout) :: ipflag + class(BndType), pointer :: packobj + integer(I4B) :: ip + + ! + ! -- Package budget summary + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_ot_bdsummary(kstp, kper, this%iout, ibudfl) + end do + + ! ! -- mover budget summary + ! if (this%inmvr > 0) then + ! call this%mvr%mvr_ot_bdsummary(ibudfl) + ! end if + + ! -- model budget summary + if (ibudfl /= 0) then + ipflag = 1 + call this%budget%budget_ot(kstp, kper, this%iout) + end if + + ! -- Write to budget csv every time step + call this%budget%writecsv(totim) + + end subroutine snf_ot_bdsummary + + !> @brief Deallocate + subroutine snf_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy + class(SnfModelType) :: this + ! -- local + integer(I4B) :: ip + class(BndType), pointer :: packobj + ! + ! -- Internal flow packages deallocate + call this%dis%dis_da() + call this%mmr%mmr_da() + call this%oc%oc_da() + call this%budget%budget_da() + ! + ! -- Internal package objects + deallocate (this%dis) + deallocate (this%budget) + deallocate (this%oc) + ! + ! -- Boundary packages + do ip = 1, this%bndlist%Count() + packobj => GetBndFromList(this%bndlist, ip) + call packobj%bnd_da() + deallocate (packobj) + end do + ! + ! -- Scalars + call mem_deallocate(this%inmmr) + call mem_deallocate(this%inoc) + ! + ! -- Arrays + ! + ! -- NumericalModelType + call this%NumericalModelType%model_da() + ! + ! -- return + return + end subroutine snf_da + + !> @brief Create boundary condition packages for this model + !< + subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, inunit, & + iout) + ! -- modules + use ConstantsModule, only: LINELENGTH + use SimModule, only: store_error + use SnfFlwModule, only: flw_create + ! -- dummy + class(SnfModelType) :: this + character(len=*), intent(in) :: filtyp + integer(I4B), intent(in) :: ipakid + integer(I4B), intent(in) :: ipaknum + character(len=*), intent(in) :: pakname + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + ! -- local + class(BndType), pointer :: packobj + class(BndType), pointer :: packobj2 + integer(I4B) :: ip + ! + ! -- This part creates the package object + select case (filtyp) + case ('FLW6') + call flw_create(packobj, ipakid, ipaknum, inunit, iout, this%name, pakname) + case default + write (errmsg, *) 'Invalid package type: ', filtyp + call store_error(errmsg, terminate=.TRUE.) + end select + ! + ! -- Check to make sure that the package name is unique, then store a + ! pointer to the package in the model bndlist + do ip = 1, this%bndlist%Count() + packobj2 => GetBndFromList(this%bndlist, ip) + if (packobj2%packName == pakname) then + write (errmsg, '(a,a)') 'Cannot create package. Package name '// & + 'already exists: ', trim(pakname) + call store_error(errmsg, terminate=.TRUE.) + end if + end do + call AddBndToList(this%bndlist, packobj) + ! + ! -- return + return + end subroutine package_create + + !> @brief Check to make sure required input files have been specified + subroutine ftype_check(this, namefile_obj, indis) + ! -- modules + use ConstantsModule, only: LINELENGTH + use SimModule, only: store_error, count_errors + use NameFileModule, only: NameFileType + ! -- dummy + class(SnfModelType) :: this + type(NameFileType), intent(in) :: namefile_obj + integer(I4B), intent(in) :: indis + ! -- local + integer(I4B) :: i, iu + character(len=LENFTYPE), dimension(4) :: nodupftype = & + (/'DISL6', 'MMR6 ', & + 'OC6 ', 'OBS6 '/) +! ------------------------------------------------------------------------------ + ! + ! -- Check for DISL, and MMR. Stop if not present. + if (indis == 0) then + write (errmsg, '(1x,a)') & + 'Discretization (DISL6) Package not specified.' + call store_error(errmsg) + end if + if (this%inmmr == 0) then + write (errmsg, '(1x,a)') & + 'Muskingum Manning Routing (MMR6) Package not specified.' + call store_error(errmsg) + end if + if (count_errors() > 0) then + write (errmsg, '(1x,a)') 'One or more required package(s) not specified.' + call store_error(errmsg) + end if + ! + ! -- Check to make sure that some SNF packages are not specified more + ! than once + do i = 1, size(nodupftype) + call namefile_obj%get_unitnumber(trim(nodupftype(i)), iu, 0) + if (iu > 0) then + write (errmsg, '(1x, a, a, a)') & + 'Duplicate entries for FTYPE ', trim(nodupftype(i)), & + ' not allowed for SNF Model.' + call store_error(errmsg) + end if + end do + ! + ! -- Stop if errors + if (count_errors() > 0) then + write (errmsg, '(a, a)') 'Error occurred while reading file: ', & + trim(namefile_obj%filename) + call store_error(errmsg, terminate=.TRUE.) + end if + ! + ! -- return + return + end subroutine ftype_check + + !> @brief Load input context for supported package + !< + subroutine snf_load_input_context(this, filtyp, modelname, pkgname, inunit, & + iout, ipaknum) + ! -- modules + use IdmMf6FileLoaderModule, only: input_load + ! -- dummy + class(SnfModelType) :: this + character(len=*), intent(in) :: filtyp + character(len=*), intent(in) :: modelname + character(len=*), intent(in) :: pkgname + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + integer(I4B), optional, intent(in) :: ipaknum + ! -- local + ! + ! -- only load if there is a file to read + if (inunit <= 0) return + ! + ! -- Load model package input to input context + select case (filtyp) + case ('MMR6') + call input_load('MMR6', 'SNF', 'MMR', modelname, pkgname, inunit, iout) + case default + call this%NumericalModelType%load_input_context(filtyp, modelname, & + pkgname, inunit, iout, & + ipaknum) + end select + ! + ! -- return + return + end subroutine snf_load_input_context + +end module SnfModule \ No newline at end of file diff --git a/src/Model/StreamNetworkFlow/snf1disl1.f90 b/src/Model/StreamNetworkFlow/snf1disl1.f90 new file mode 100644 index 00000000000..06a81309d95 --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1disl1.f90 @@ -0,0 +1,1059 @@ +module SnfDislModule + + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: LENMEMPATH, LENVARNAME, DZERO, DONE, LINELENGTH + use SimVariablesModule, only: errmsg, warnmsg + use MemoryHelperModule, only: create_mem_path + use MemoryManagerModule, only: mem_allocate + use SimModule, only: count_errors, store_error, store_error_unit, & + store_warning + use BaseDisModule, only: DisBaseType + use DislGeom, only: calcdist, partialdist + + implicit none + + private + public :: disl_cr + public :: SnfDislType + + type, extends(DisBaseType) :: SnfDislType + integer(I4B), pointer :: nvert => null() !< number of x,y vertices + real(DP), pointer :: convlength => null() !< conversion factor for length + real(DP), pointer :: convtime => null() !< conversion factor for time + real(DP), dimension(:), pointer, contiguous :: segment_length => null() !< length of each segment + integer(I4B), dimension(:), pointer, contiguous :: tosegment => null() !< downstream segment index (nodes) + integer(I4B), dimension(:), pointer, contiguous :: idomain => null() !< idomain (nodes) + real(DP), dimension(:, :), pointer, contiguous :: vertices => null() !< cell vertices stored as 2d array with columns of x, y, and z + real(DP), dimension(:, :), pointer, contiguous :: cellxyz => null() !< segment midpoints stored as 2d array with columns of x, y, and z + real(DP), dimension(:), pointer, contiguous :: fdc => null() !< fdc stored as array + integer(I4B), dimension(:), pointer, contiguous :: iavert => null() !< cell vertex pointer ia array + integer(I4B), dimension(:), pointer, contiguous :: javert => null() !< cell vertex pointer ja array + contains + procedure :: disl_load + procedure :: dis_da => disl_da + procedure, public :: record_array + procedure, public :: record_srcdst_list_header + ! -- private + procedure :: allocate_scalars + procedure :: allocate_arrays + procedure :: source_options + procedure :: source_dimensions + procedure :: source_griddata + procedure :: source_vertices + procedure :: source_cell2d + procedure :: log_options + procedure :: log_dimensions + procedure :: log_griddata + procedure :: define_cellverts + procedure :: grid_finalize + !procedure :: connect + procedure :: create_connections + procedure :: write_grb + procedure :: get_nodenumber_idx1 + procedure :: nodeu_to_string + + end type SnfDislType + +contains + + !> @brief Create new disl package + !< + subroutine disl_cr(dis, name_model, inunit, iout) + use IdmMf6FileLoaderModule, only: input_load + class(DisBaseType), pointer :: dis + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in) :: iout + type(SnfDislType), pointer :: disnew + character(len=*), parameter :: fmtheader = & + "(1X, /1X, 'DISL -- LINE NETWORK DISCRETIZATION PACKAGE,', & + &' VERSION 1 : 1/20/2023 - INPUT READ FROM UNIT ', I0, //)" + allocate (disnew) + dis => disnew + call disnew%allocate_scalars(name_model) + dis%inunit = inunit + dis%iout = iout + ! + ! -- if reading from file + if (inunit > 0) then + ! + ! -- Identify package + if (iout > 0) then + write (iout, fmtheader) inunit + end if + ! + ! -- load disl + call disnew%disl_load() + end if + ! + ! -- Return + return + end subroutine disl_cr + + !> @brief Allocate scalar variables + !< + subroutine allocate_scalars(this, name_model) + ! -- modules + use MemoryManagerModule, only: mem_allocate + use ConstantsModule, only: DONE + ! -- dummy + class(SnfDislType) :: this + character(len=*), intent(in) :: name_model + ! + ! -- Allocate parent scalars + call this%DisBaseType%allocate_scalars(name_model) + ! + ! -- Allocate + call mem_allocate(this%nvert, 'NVERT', this%memoryPath) + call mem_allocate(this%convlength, 'CONVLENGTH', this%memoryPath) + call mem_allocate(this%convtime, 'CONVTIME', this%memoryPath) + ! + ! -- Initialize + this%nvert = 0 + this%ndim = 1 + this%convlength = DONE + this%convtime = DONE + ! + ! -- Return + return + end subroutine allocate_scalars + + subroutine disl_load(this) + ! -- dummy + class(SnfDislType) :: this + ! -- locals + ! + ! -- source input data + call this%source_options() + call this%source_dimensions() + call this%source_griddata() + + ! create connectivity using tosegment + call this%create_connections() + + ! If vertices provided by user, read and store vertices + if (this%nvert > 0) then + call this%source_vertices() + call this%source_cell2d() + end if + + ! finalize the grid + call this%grid_finalize() + ! + ! -- Return + return + end subroutine disl_load + + !> @brief Copy options from IDM into package + !< + subroutine source_options(this) + ! -- modules + use KindModule, only: LGP + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + use SnfDislInputModule, only: SnfDislParamFoundType + ! -- dummy + class(SnfDislType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + character(len=LENVARNAME), dimension(3) :: lenunits = & + &[character(len=LENVARNAME) :: 'FEET', 'METERS', 'CENTIMETERS'] + type(SnfDislParamFoundType) :: found + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'DISL', idm_context) + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%lenuni, 'LENGTH_UNITS', idmMemoryPath, lenunits, & + found%length_units) + call mem_set_value(this%convlength, 'CONVLENGTH', idmMemoryPath, found%length_convert) + call mem_set_value(this%convtime, 'CONVTIME', idmMemoryPath, found%time_convert) + call mem_set_value(this%nogrb, 'NOGRB', idmMemoryPath, found%nogrb) + call mem_set_value(this%xorigin, 'XORIGIN', idmMemoryPath, found%xorigin) + call mem_set_value(this%yorigin, 'YORIGIN', idmMemoryPath, found%yorigin) + call mem_set_value(this%angrot, 'ANGROT', idmMemoryPath, found%angrot) + ! + ! -- log values to list file + if (this%iout > 0) then + call this%log_options(found) + end if + ! + ! -- Return + return + end subroutine source_options + + !> @brief Write user options to list file + !< + subroutine log_options(this, found) + use SnfDislInputModule, only: SnfDislParamFoundType + class(SnfDislType) :: this + type(SnfDislParamFoundType), intent(in) :: found + + write (this%iout, '(1x,a)') 'Setting Discretization Options' + + if (found%length_units) then + write (this%iout, '(4x,a,i0)') 'Model length unit [0=UND, 1=FEET, & + &2=METERS, 3=CENTIMETERS] set as ', this%lenuni + end if + + if (found%nogrb) then + write (this%iout, '(4x,a,i0)') 'Binary grid file [0=GRB, 1=NOGRB] & + &set as ', this%nogrb + end if + + if (found%xorigin) then + write (this%iout, '(4x,a,G0)') 'XORIGIN = ', this%xorigin + end if + + if (found%yorigin) then + write (this%iout, '(4x,a,G0)') 'YORIGIN = ', this%yorigin + end if + + if (found%angrot) then + write (this%iout, '(4x,a,G0)') 'ANGROT = ', this%angrot + end if + + if (found%length_convert) then + write (this%iout, '(4x,a,G0)') 'LENGTH_CONVERSION = ', this%convlength + end if + + if (found%time_convert) then + write (this%iout, '(4x,a,G0)') 'TIME_CONVERSION = ', this%convtime + end if + + write (this%iout, '(1x,a,/)') 'End Setting Discretization Options' + + end subroutine log_options + + !> @brief Copy dimensions from IDM into package + !< + subroutine source_dimensions(this) + use KindModule, only: LGP + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + use SnfDislInputModule, only: SnfDislParamFoundType + ! -- dummy + class(SnfDislType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + integer(I4B) :: n + type(SnfDislParamFoundType) :: found + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'DISL', idm_context) + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%nodes, 'NODES', idmMemoryPath, found%nodes) + call mem_set_value(this%nvert, 'NVERT', idmMemoryPath, found%nvert) + ! + ! -- for now assume nodes = nodesuser + this%nodesuser = this%nodes + ! + ! -- log simulation values + if (this%iout > 0) then + call this%log_dimensions(found) + end if + ! + ! -- verify dimensions were set + if (this%nodesuser < 1) then + call store_error( & + 'NODES was not specified or was specified incorrectly.') + call store_error_unit(this%inunit) + end if + if (this%nvert < 1) then + call store_warning( & + 'NVERT was not specified or was specified as zero. The & + &VERTICES and CELL2D blocks will not be read for the DISL6 & + &Package in model ' // trim(this%memoryPath) // '.') + end if + ! + ! -- Allocate non-reduced vectors for disl + call mem_allocate(this%segment_length, this%nodesuser, 'SEGMENT_LENGTH', this%memoryPath) + call mem_allocate(this%tosegment, this%nodesuser, 'TOSEGMENT', this%memoryPath) + call mem_allocate(this%idomain, this%nodesuser, 'IDOMAIN', this%memoryPath) + ! + ! -- Allocate vertices array + if (this%nvert > 0) then + call mem_allocate(this%vertices, 3, this%nvert, 'VERTICES', this%memoryPath) + call mem_allocate(this%fdc, this%nodesuser, 'FDC', this%memoryPath) + call mem_allocate(this%cellxyz, 3, this%nodesuser, 'CELLXYZ', this%memoryPath) + end if + ! + ! -- initialize all cells to be active (idomain = 1) + do n = 1, this%nodesuser + this%segment_length(n) = DZERO + this%tosegment(n) = 0 + this%idomain(n) = 1 + end do + ! + ! -- Return + return + end subroutine source_dimensions + + !> @brief Write dimensions to list file + !< + subroutine log_dimensions(this, found) + use SnfDislInputModule, only: SnfDislParamFoundType + class(SnfDislType) :: this + type(SnfDislParamFoundType), intent(in) :: found + + write (this%iout, '(1x,a)') 'Setting Discretization Dimensions' + + if (found%nodes) then + write (this%iout, '(4x,a,i0)') 'NODES = ', this%nodesuser + end if + + if (found%nvert) then + write (this%iout, '(4x,a,i0)') 'NVERT = ', this%nvert + end if + + write (this%iout, '(1x,a,/)') 'End Setting Discretization Dimensions' + + end subroutine log_dimensions + + subroutine source_griddata(this) + ! -- modules + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + use SnfDislInputModule, only: SnfDislParamFoundType + ! -- dummy + class(SnfDislType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + type(SnfDislParamFoundType) :: found + ! -- formats + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'DISL', idm_context) + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%segment_length, 'SEGMENT_LENGTH', idmMemoryPath, & + found%segment_length) + call mem_set_value(this%tosegment, 'TOSEGMENT', idmMemoryPath, & + found%tosegment) + call mem_set_value(this%idomain, 'IDOMAIN', idmMemoryPath, found%idomain) + ! + ! -- log simulation values + if (this%iout > 0) then + call this%log_griddata(found) + end if + ! + ! -- Return + return + end subroutine source_griddata + + !> @brief Write griddata found to list file + !< + subroutine log_griddata(this, found) + use SnfDislInputModule, only: SnfDislParamFoundType + class(SnfDislType) :: this + type(SnfDislParamFoundType), intent(in) :: found + + write (this%iout, '(1x,a)') 'Setting Discretization Griddata' + + if (found%segment_length) then + write (this%iout, '(4x,a)') 'SEGMENT_LENGTH set from input file' + end if + + if (found%tosegment) then + write (this%iout, '(4x,a)') 'TOSEGMENT set from input file' + end if + + if (found%idomain) then + write (this%iout, '(4x,a)') 'IDOMAIN set from input file' + end if + + write (this%iout, '(1x,a,/)') 'End Setting Discretization Griddata' + + end subroutine log_griddata + + !> @brief Copy vertex information from input data context + !! to model context + !< + subroutine source_vertices(this) + ! -- modules + use MemoryManagerModule, only: mem_setptr + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + ! -- dummy + class(SnfDislType) :: this + ! -- local + integer(I4B) :: i + character(len=LENMEMPATH) :: idmMemoryPath + real(DP), dimension(:), contiguous, pointer :: vert_x => null() + real(DP), dimension(:), contiguous, pointer :: vert_y => null() + real(DP), dimension(:), contiguous, pointer :: vert_z => null() + ! -- formats +! ------------------------------------------------------------------------------ + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'DISL', idm_context) + ! + ! -- set pointers to memory manager input arrays + call mem_setptr(vert_x, 'XV', idmMemoryPath) + call mem_setptr(vert_y, 'YV', idmMemoryPath) + call mem_setptr(vert_z, 'ZV', idmMemoryPath) + ! + ! -- set vertices 3d array + if (associated(vert_x) .and. associated(vert_y) .and. & + associated(vert_z)) then + do i = 1, this%nvert + this%vertices(1, i) = vert_x(i) + this%vertices(2, i) = vert_y(i) + this%vertices(3, i) = vert_z(i) + end do + else + call store_error('Required Vertex arrays not found.') + end if + ! + ! -- log + if (this%iout > 0) then + write (this%iout, '(1x,a)') 'Setting Discretization Vertices' + write (this%iout, '(1x,a,/)') 'End setting discretization vertices' + end if + ! + ! -- Return + return + end subroutine source_vertices + + !> @brief Copy cell2d information from input data context + !! to model context + !< + subroutine source_cell2d(this) + ! -- modules + use MemoryHelperModule, only: create_mem_path + use MemoryManagerModule, only: mem_setptr + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + ! -- dummy + class(SnfDislType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + integer(I4B), dimension(:), contiguous, pointer :: icell2d => null() + integer(I4B), dimension(:), contiguous, pointer :: ncvert => null() + integer(I4B), dimension(:), contiguous, pointer :: icvert => null() + real(DP), dimension(:), contiguous, pointer :: fdc => null() + integer(I4B) :: i + ! -- formats + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'DISL', idm_context) + ! + ! -- set pointers to input path ncvert and icvert + call mem_setptr(icell2d, 'ICELL2D', idmMemoryPath) + call mem_setptr(ncvert, 'NCVERT', idmMemoryPath) + call mem_setptr(icvert, 'ICVERT', idmMemoryPath) + ! + ! -- + if (associated(icell2d) .and. associated(ncvert) & + .and. associated(icvert)) then + call this%define_cellverts(icell2d, ncvert, icvert) + else + call store_error('Required cell vertex arrays not found.') + end if + ! + ! -- set pointers to cell center arrays + call mem_setptr(fdc, 'FDC', idmMemoryPath) + ! + ! -- set fractional distance to cell center + if (associated(fdc)) then + do i = 1, this%nodesuser + this%fdc(i) = fdc(i) + end do + call calculate_cellxyz(this%vertices, this%fdc, this%iavert, & + this%javert, this%cellxyz) + else + call store_error('Required fdc array not found.') + end if + ! + ! -- log + if (this%iout > 0) then + write (this%iout, '(1x,a)') 'Setting Discretization CELL2D' + write (this%iout, '(1x,a,/)') 'End Setting Discretization CELL2D' + end if + ! + ! -- Return + return + end subroutine source_cell2d + + !> @brief Construct the iavert and javert integer vectors which + !! are compressed sparse row index arrays that relate the vertices + !! to segments + !< + subroutine define_cellverts(this, icell2d, ncvert, icvert) + ! -- modules + use SparseModule, only: sparsematrix + ! -- dummy + class(SnfDislType) :: this + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icell2d + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: ncvert + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: icvert + ! -- locals + type(sparsematrix) :: vert_spm + integer(I4B) :: i, j, ierr + integer(I4B) :: icv_idx, startvert, maxnnz = 2 +! ------------------------------------------------------------------------------ + ! + ! -- initialize sparse matrix + call vert_spm%init(this%nodesuser, this%nvert, maxnnz) + ! + ! -- add sparse matrix connections from input memory paths + icv_idx = 1 + do i = 1, this%nodesuser + if (icell2d(i) /= i) call store_error('ICELL2D input sequence violation.') + do j = 1, ncvert(i) + call vert_spm%addconnection(i, icvert(icv_idx), 0) + if (j == 1) then + startvert = icvert(icv_idx) + end if + icv_idx = icv_idx + 1 + end do + end do + ! + ! -- allocate and fill iavert and javert + call mem_allocate(this%iavert, this%nodesuser + 1, 'IAVERT', this%memoryPath) + call mem_allocate(this%javert, vert_spm%nnz, 'JAVERT', this%memoryPath) + call vert_spm%filliaja(this%iavert, this%javert, ierr) + call vert_spm%destroy() + ! + ! -- Return + return + end subroutine define_cellverts + + !> @brief Calculate x, y, z coordinates of segment midpoint + !< + subroutine calculate_cellxyz(vertices, fdc, iavert, javert, cellxyz) + ! -- dummy + real(DP), dimension(:, :), intent(in) :: vertices !< 2d array of vertices with x, y, and z as columns + real(DP), dimension(:), intent(in) :: fdc !< fractional distance to segment midpoint (normally 0.5) + integer(I4B), dimension(:), intent(in) :: iavert !< csr mapping of vertices to cell segments + integer(I4B), dimension(:), intent(in) :: javert !< csr mapping of vertices to cell segments + real(DP), dimension(:, :), intent(inout) :: cellxyz !< 2d array of segment midpoint with x, y, and z as columns + ! -- local + integer(I4B) :: nodes !< number of nodes + integer(I4B) :: n !< node index + integer(I4B) :: j !< vertex index + integer(I4B) :: iv0 !< index for line segment start + integer(I4B) :: iv1 !< index for linen segment end + integer(I4B) :: ixyz !< x, y, z column index + real(DP) :: segment_length !< segment length = sum of individual line segments + real(DP) :: fd0 !< fractional distance to start of this line segment + real(DP) :: fd1 !< fractional distance to end fo this line segment + real(DP) :: fd !< fractional distance where midpoint (defined by fdc) is located + real(DP) :: d !< distance + + nodes = size(iavert) - 1 + do n = 1, nodes + + ! calculate length of this segment + segment_length = DZERO + do j = iavert(n), iavert(n + 1) - 2 + segment_length = segment_length + & + calcdist(vertices, javert(j), javert(j + 1)) + end do + + ! find vertices that span midpoint + iv0 = 0 + iv1 = 0 + fd0 = DZERO + do j = iavert(n), iavert(n + 1) - 2 + d = calcdist(vertices, javert(j), javert(j + 1)) + fd1 = fd0 + d / segment_length + + ! if true, then we know the midpoint is some fractional distance (fd) + ! from vertex j to vertex j + 1 + if (fd1 >= fdc(n)) then + iv0 = javert(j) + iv1 = javert(j + 1) + fd = (fdc(n) - fd0) / (fd1 - fd0) + exit + end if + fd0 = fd1 + end do + + ! find x, y, z position of point on line + do ixyz = 1, 3 + cellxyz(ixyz, n) = (DONE - fd) * vertices(ixyz, iv0) + & + fd * vertices(ixyz, iv1) + end do + + end do + end subroutine calculate_cellxyz + + !> @brief Finalize grid construction + !< + subroutine grid_finalize(this) + ! -- modules + use SimModule, only: ustop, count_errors, store_error + use ConstantsModule, only: LINELENGTH, DZERO, DONE + ! -- dummy + class(SnfDislType) :: this + ! -- locals + integer(I4B) :: node, noder, k + ! -- formats + character(len=*), parameter :: fmtdz = & + "('ERROR. CELL (',i0,',',i0,') THICKNESS <= 0. ', " // & + "'TOP, BOT: ',2(1pg24.15))" + character(len=*), parameter :: fmtnr = & + "(/1x, 'THE SPECIFIED IDOMAIN RESULTS IN A REDUCED NUMBER OF CELLS.'," // & + "/1x, 'NUMBER OF USER NODES: ',I7," // & + "/1X, 'NUMBER OF NODES IN SOLUTION: ', I7, //)" + ! -- data + ! + ! -- count active cells + this%nodes = 0 + do k = 1, this%nodesuser + if(this%idomain(k) > 0) this%nodes = this%nodes + 1 + enddo + ! + ! -- Check to make sure nodes is a valid number + if (this%nodes == 0) then + call store_error('ERROR. MODEL DOES NOT HAVE ANY ACTIVE NODES.') + call store_error('MAKE SURE IDOMAIN ARRAY HAS SOME VALUES GREATER & + &THAN ZERO.') + call this%parser%StoreErrorUnit() + call ustop() + end if + + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + call ustop() + endif + ! + ! -- Array size is now known, so allocate + call this%allocate_arrays() + ! + ! -- Fill the nodereduced array with the reduced nodenumber, or + ! a negative number to indicate it is a pass-through cell, or + ! a zero to indicate that the cell is excluded from the + ! solution. + if(this%nodes < this%nodesuser) then + node = 1 + noder = 1 + do k = 1, this%nodesuser + if(this%idomain(k) > 0) then + this%nodereduced(node) = noder + noder = noder + 1 + elseif(this%idomain(k) < 0) then + this%nodereduced(node) = -1 + else + this%nodereduced(node) = 0 + endif + node = node + 1 + enddo + endif + ! + ! -- allocate and fill nodeuser if a reduced grid + if(this%nodes < this%nodesuser) then + node = 1 + noder = 1 + do k = 1, this%nodesuser + if(this%idomain(k) > 0) then + this%nodeuser(noder) = node + noder = noder + 1 + endif + node = node + 1 + enddo + endif + + ! -- Return + return + end subroutine grid_finalize + + subroutine allocate_arrays(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy + class(SnfDislType) :: this + ! + ! -- Allocate arrays in DisBaseType (mshape, top, bot, area) + call this%DisBaseType%allocate_arrays() + ! + ! -- Allocate arrays + if(this%nodes < this%nodesuser) then + call mem_allocate(this%nodeuser, this%nodes, 'NODEUSER', this%memoryPath) + call mem_allocate(this%nodereduced, this%nodesuser, 'NODEREDUCED', & + this%memoryPath) + else + call mem_allocate(this%nodeuser, 1, 'NODEUSER', this%memoryPath) + call mem_allocate(this%nodereduced, 1, 'NODEREDUCED', this%memoryPath) + endif + ! + ! -- Initialize + this%mshape(1) = this%nodesuser + ! + ! -- Return + return + end subroutine allocate_arrays + + subroutine create_connections(this) + ! -- modules + ! -- dummy + class(SnfDislType) :: this + ! -- local + integer(I4B) :: nrsize + ! + ! -- create and fill the connections object + nrsize = 0 + if(this%nodes < this%nodesuser) nrsize = this%nodes + ! + ! -- Allocate connections object + allocate(this%con) + + call this%con%dislconnections(this%name_model, this%tosegment) + this%nja = this%con%nja + this%njas = this%con%njas + ! + ! + ! -- return + return + end subroutine create_connections + + !> @brief Write binary grid file + !< + subroutine write_grb(this, icelltype) + ! -- modules + use InputOutputModule, only: getunit, openfile + use OpenSpecModule, only: access, form + ! -- dummy + class(SnfDislType) :: this + integer(I4B), dimension(:), intent(in) :: icelltype + ! -- local + integer(I4B) :: i, iunit, ntxt + integer(I4B), parameter :: lentxt = 100 + character(len=50) :: txthdr + character(len=lentxt) :: txt + character(len=LINELENGTH) :: fname + character(len=*), parameter :: fmtgrdsave = & + "(4X,'BINARY GRID INFORMATION WILL BE WRITTEN TO:', & + &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" + ! + ! -- Initialize + ntxt = 9 + if (this%nvert > 0) ntxt = ntxt + 6 + ! + ! -- Open the file + inquire (unit=this%inunit, name=fname) + fname = trim(fname)//'.grb' + iunit = getunit() + write (this%iout, fmtgrdsave) iunit, trim(adjustl(fname)) + call openfile(iunit, this%iout, trim(adjustl(fname)), 'DATA(BINARY)', & + form, access, 'REPLACE') + ! + ! -- write header information + write (txthdr, '(a)') 'GRID DISL' + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a)') 'VERSION 1' + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a, i0)') 'NTXT ', ntxt + txthdr(50:50) = new_line('a') + write (iunit) txthdr + write (txthdr, '(a, i0)') 'LENTXT ', lentxt + txthdr(50:50) = new_line('a') + write (iunit) txthdr + ! + ! -- write variable definitions + write (txt, '(3a, i0)') 'NODES ', 'INTEGER ', 'NDIM 0 # ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'NJA ', 'INTEGER ', 'NDIM 0 # ', this%con%nja + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg24.15)') 'XORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%xorigin + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg24.15)') 'YORIGIN ', 'DOUBLE ', 'NDIM 0 # ', this%yorigin + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, 1pg24.15)') 'ANGROT ', 'DOUBLE ', 'NDIM 0 # ', this%angrot + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IA ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1 + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', this%con%nja + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + ! + ! -- if vertices have been read then write additional header information + if (this%nvert > 0) then + write (txt, '(3a, i0)') 'VERTICES ', 'DOUBLE ', 'NDIM 2 3 ', this%nvert + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'CELLX ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'CELLY ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'CELLZ ', 'DOUBLE ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'IAVERT ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + 1 + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + write (txt, '(3a, i0)') 'JAVERT ', 'INTEGER ', 'NDIM 1 ', size(this%javert) + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt + end if + ! + ! -- write data + write (iunit) this%nodesuser ! nodes + write (iunit) this%nja ! nja + write (iunit) this%xorigin ! xorigin + write (iunit) this%yorigin ! yorigin + write (iunit) this%angrot ! angrot + write (iunit) this%con%iausr ! ia + write (iunit) this%con%jausr ! ja + write (iunit) icelltype ! icelltype + write (iunit) this%idomain ! idomain + ! + ! -- if vertices have been read then write additional data + if (this%nvert > 0) then + write (iunit) this%vertices ! vertices + write (iunit) (this%cellxyz(1, i), i=1, this%nodesuser) ! cellx + write (iunit) (this%cellxyz(2, i), i=1, this%nodesuser) ! celly + write (iunit) (this%cellxyz(3, i), i=1, this%nodesuser) ! cellz + write (iunit) this%iavert ! iavert + write (iunit) this%javert ! javert + end if + ! + ! -- Close the file + close (iunit) + ! + ! -- return + return + end subroutine write_grb + + !> + !! Return a nodenumber from the user specified node number with an + !! option to perform a check. This subroutine can be overridden by + !! child classes to perform mapping to a model node number + !< + function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) + class(SnfDislType), intent(in) :: this + integer(I4B), intent(in) :: nodeu + integer(I4B), intent(in) :: icheck + integer(I4B) :: nodenumber + ! + if (icheck /= 0) then + if (nodeu < 1 .or. nodeu > this%nodes) then + write (errmsg, '(a,i10)') & + 'Nodenumber less than 1 or greater than nodes:', nodeu + call store_error(errmsg) + end if + end if + ! + ! -- set node number based on wheter it is reduced or not + if (this%nodes == this%nodesuser) then + nodenumber = nodeu + else + nodenumber = this%nodereduced(nodeu) + end if + ! + ! -- return + return + end function get_nodenumber_idx1 + + subroutine nodeu_to_string(this, nodeu, str) + ! -- dummy + class(SnfDislType) :: this + integer(I4B), intent(in) :: nodeu + character(len=*), intent(inout) :: str + ! -- local + character(len=10) :: nstr + ! + write (nstr, '(i0)') nodeu + str = '('//trim(adjustl(nstr))//')' + ! + ! -- return + return + end subroutine nodeu_to_string + + subroutine disl_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + use MemoryManagerExtModule, only: memorylist_remove + use SimVariablesModule, only: idm_context + ! -- dummy + class(SnfDislType) :: this + ! -- local + logical(LGP) :: deallocate_vertices + ! + ! -- Deallocate idm memory + call memorylist_remove(this%name_model, 'DISL', idm_context) + call memorylist_remove(component=this%name_model, & + context=idm_context) + ! + ! -- scalars + deallocate_vertices = (this%nvert > 0) + call mem_deallocate(this%nvert) + call mem_deallocate(this%convlength) + call mem_deallocate(this%convtime) + ! + ! -- arrays + call mem_deallocate(this%nodeuser) + call mem_deallocate(this%nodereduced) + call mem_deallocate(this%segment_length) + call mem_deallocate(this%tosegment) + call mem_deallocate(this%idomain) + ! + ! -- cdl hack for arrays for vertices and cell2d blocks + if (deallocate_vertices) then + call mem_deallocate(this%vertices) + call mem_deallocate(this%fdc) + call mem_deallocate(this%cellxyz) + call mem_deallocate(this%iavert) + call mem_deallocate(this%javert) + end if + ! + ! -- DisBaseType deallocate + call this%DisBaseType%dis_da() + ! + ! -- Return + return + end subroutine disl_da + + + !> @brief Record a double precision array + !! + !! Record a double precision array. The array will be + !! printed to an external file and/or written to an unformatted external file + !! depending on the argument specifications. + !< + subroutine record_array(this, darray, iout, iprint, idataun, aname, & + cdatafmp, nvaluesp, nwidthp, editdesc, dinact) + ! -- modules + use TdisModule, only: kstp, kper, pertim, totim, delt + use InputOutputModule, only: ulasav, ulaprufw, ubdsv1 + ! -- dummy + class(SnfDislType), intent(inout) :: this + real(DP), dimension(:), pointer, contiguous, intent(inout) :: darray !< double precision array to record + integer(I4B), intent(in) :: iout !< unit number for ascii output + integer(I4B), intent(in) :: iprint !< flag indicating whether or not to print the array + integer(I4B), intent(in) :: idataun !< unit number to which the array will be written in binary + character(len=*), intent(in) :: aname !< text descriptor of the array + character(len=*), intent(in) :: cdatafmp ! fortran format for writing the array + integer(I4B), intent(in) :: nvaluesp !< number of values per line for printing + integer(I4B), intent(in) :: nwidthp !< width of the number for printing + character(len=*), intent(in) :: editdesc !< format type (I, G, F, S, E) + real(DP), intent(in) :: dinact !< double precision value to use for cells that are excluded from model domain + ! -- local + integer(I4B) :: k, ifirst + integer(I4B) :: nlay + integer(I4B) :: nrow + integer(I4B) :: ncol + integer(I4B) :: nval + integer(I4B) :: nodeu, noder + integer(I4B) :: istart, istop + real(DP), dimension(:), pointer, contiguous :: dtemp + ! -- formats + character(len=*), parameter :: fmthsv = & + "(1X,/1X,a,' WILL BE SAVED ON UNIT ',I4, & + &' AT END OF TIME STEP',I5,', STRESS PERIOD ',I4)" + ! + ! -- set variables + nlay = 1 + nrow = 1 + ncol = this%mshape(1) + ! + ! -- If this is a reduced model, then copy the values from darray into + ! dtemp. + if (this%nodes < this%nodesuser) then + nval = this%nodes + dtemp => this%dbuff + do nodeu = 1, this%nodesuser + noder = this%get_nodenumber(nodeu, 0) + if (noder <= 0) then + dtemp(nodeu) = dinact + cycle + end if + dtemp(nodeu) = darray(noder) + end do + else + nval = this%nodes + dtemp => darray + end if + ! + ! -- Print to iout if iprint /= 0 + if (iprint /= 0) then + istart = 1 + do k = 1, nlay + istop = istart + nrow * ncol - 1 + call ulaprufw(ncol, nrow, kstp, kper, k, iout, dtemp(istart:istop), & + aname, cdatafmp, nvaluesp, nwidthp, editdesc) + istart = istop + 1 + end do + end if + ! + ! -- Save array to an external file. + if (idataun > 0) then + ! -- write to binary file by layer + ifirst = 1 + istart = 1 + do k = 1, nlay + istop = istart + nrow * ncol - 1 + if (ifirst == 1) write (iout, fmthsv) & + trim(adjustl(aname)), idataun, & + kstp, kper + ifirst = 0 + call ulasav(dtemp(istart:istop), aname, kstp, kper, & + pertim, totim, ncol, nrow, k, idataun) + istart = istop + 1 + end do + elseif (idataun < 0) then + ! + ! -- write entire array as one record + call ubdsv1(kstp, kper, aname, -idataun, dtemp, ncol, nrow, nlay, & + iout, delt, pertim, totim) + end if + ! + ! -- return + return + end subroutine record_array + + !> @brief Record list header using ubdsv06 + !< + subroutine record_srcdst_list_header(this, text, textmodel, textpackage, & + dstmodel, dstpackage, naux, auxtxt, & + ibdchn, nlist, iout) + ! -- module + use TdisModule, only: kstp, kper, pertim, totim, delt + use InputOutputModule, only: ubdsv06 + ! -- dummy + class(SnfDislType) :: this + character(len=16), intent(in) :: text + character(len=16), intent(in) :: textmodel + character(len=16), intent(in) :: textpackage + character(len=16), intent(in) :: dstmodel + character(len=16), intent(in) :: dstpackage + integer(I4B), intent(in) :: naux + character(len=16), dimension(:), intent(in) :: auxtxt + integer(I4B), intent(in) :: ibdchn + integer(I4B), intent(in) :: nlist + integer(I4B), intent(in) :: iout + ! -- local + integer(I4B) :: nlay, nrow, ncol + ! + nlay = 1 + nrow = 1 + ncol = this%mshape(1) + ! + ! -- Use ubdsv06 to write list header + call ubdsv06(kstp, kper, text, textmodel, textpackage, dstmodel, dstpackage, & + ibdchn, naux, auxtxt, ncol, nrow, nlay, & + nlist, iout, delt, pertim, totim) + ! + ! -- return + return + end subroutine record_srcdst_list_header + +end module SnfDislModule \ No newline at end of file diff --git a/src/Model/StreamNetworkFlow/snf1disl1idm.f90 b/src/Model/StreamNetworkFlow/snf1disl1idm.f90 new file mode 100644 index 00000000000..7d03f5124e4 --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1disl1idm.f90 @@ -0,0 +1,447 @@ +module SnfDislInputModule + use InputDefinitionModule, only: InputParamDefinitionType, & + InputBlockDefinitionType + private + public snf_disl_param_definitions + public snf_disl_aggregate_definitions + public snf_disl_block_definitions + public SnfDislParamFoundType + + type SnfDislParamFoundType + logical :: length_units = .false. + logical :: length_convert = .false. + logical :: time_convert = .false. + logical :: nogrb = .false. + logical :: xorigin = .false. + logical :: yorigin = .false. + logical :: angrot = .false. + logical :: nodes = .false. + logical :: nvert = .false. + logical :: segment_length = .false. + logical :: tosegment = .false. + logical :: idomain = .false. + logical :: iv = .false. + logical :: xv = .false. + logical :: yv = .false. + logical :: zv = .false. + logical :: icell2d = .false. + logical :: fdc = .false. + logical :: ncvert = .false. + logical :: icvert = .false. + end type SnfDislParamFoundType + + type(InputParamDefinitionType), parameter :: & + snfdisl_length_units = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'LENGTH_UNITS', & ! tag name + 'LENGTH_UNITS', & ! fortran variable + 'STRING', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_length_convert = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'LENGTH_CONVERT', & ! tag name + 'LENGTH_CONVERT', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_time_convert = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'TIME_CONVERT', & ! tag name + 'TIME_CONVERT', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_nogrb = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'NOGRB', & ! tag name + 'NOGRB', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_xorigin = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'XORIGIN', & ! tag name + 'XORIGIN', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_yorigin = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'YORIGIN', & ! tag name + 'YORIGIN', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_angrot = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'OPTIONS', & ! block + 'ANGROT', & ! tag name + 'ANGROT', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_nodes = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'DIMENSIONS', & ! block + 'NODES', & ! tag name + 'NODES', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_nvert = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'DIMENSIONS', & ! block + 'NVERT', & ! tag name + 'NVERT', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_segment_length = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'GRIDDATA', & ! block + 'SEGMENT_LENGTH', & ! tag name + 'SEGMENT_LENGTH', & ! fortran variable + 'DOUBLE1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_tosegment = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'GRIDDATA', & ! block + 'TOSEGMENT', & ! tag name + 'TOSEGMENT', & ! fortran variable + 'INTEGER1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_idomain = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'GRIDDATA', & ! block + 'IDOMAIN', & ! tag name + 'IDOMAIN', & ! fortran variable + 'INTEGER1D', & ! type + 'NODES', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_iv = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'VERTICES', & ! block + 'IV', & ! tag name + 'IV', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_xv = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'VERTICES', & ! block + 'XV', & ! tag name + 'XV', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_yv = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'VERTICES', & ! block + 'YV', & ! tag name + 'YV', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_zv = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'VERTICES', & ! block + 'ZV', & ! tag name + 'ZV', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_icell2d = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'CELL2D', & ! block + 'ICELL2D', & ! tag name + 'ICELL2D', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_fdc = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'CELL2D', & ! block + 'FDC', & ! tag name + 'FDC', & ! fortran variable + 'DOUBLE', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_ncvert = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'CELL2D', & ! block + 'NCVERT', & ! tag name + 'NCVERT', & ! fortran variable + 'INTEGER', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_icvert = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'CELL2D', & ! block + 'ICVERT', & ! tag name + 'ICVERT', & ! fortran variable + 'INTEGER1D', & ! type + 'NCVERT', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snf_disl_param_definitions(*) = & + [ & + snfdisl_length_units, & + snfdisl_length_convert, & + snfdisl_time_convert, & + snfdisl_nogrb, & + snfdisl_xorigin, & + snfdisl_yorigin, & + snfdisl_angrot, & + snfdisl_nodes, & + snfdisl_nvert, & + snfdisl_segment_length, & + snfdisl_tosegment, & + snfdisl_idomain, & + snfdisl_iv, & + snfdisl_xv, & + snfdisl_yv, & + snfdisl_zv, & + snfdisl_icell2d, & + snfdisl_fdc, & + snfdisl_ncvert, & + snfdisl_icvert & + ] + + type(InputParamDefinitionType), parameter :: & + snfdisl_vertices = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'VERTICES', & ! block + 'VERTICES', & ! tag name + 'VERTICES', & ! fortran variable + 'RECARRAY IV XV YV ZV', & ! type + 'NVERT', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfdisl_cell2d = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'DISL', & ! subcomponent + 'CELL2D', & ! block + 'CELL2D', & ! tag name + 'CELL2D', & ! fortran variable + 'RECARRAY ICELL2D FDC NCVERT ICVERT', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snf_disl_aggregate_definitions(*) = & + [ & + snfdisl_vertices, & + snfdisl_cell2d & + ] + + type(InputBlockDefinitionType), parameter :: & + snf_disl_block_definitions(*) = & + [ & + InputBlockDefinitionType( & + 'OPTIONS', & ! blockname + .false., & ! required + .false. & ! aggregate + ), & + InputBlockDefinitionType( & + 'DIMENSIONS', & ! blockname + .true., & ! required + .false. & ! aggregate + ), & + InputBlockDefinitionType( & + 'GRIDDATA', & ! blockname + .true., & ! required + .false. & ! aggregate + ), & + InputBlockDefinitionType( & + 'VERTICES', & ! blockname + .true., & ! required + .true. & ! aggregate + ), & + InputBlockDefinitionType( & + 'CELL2D', & ! blockname + .true., & ! required + .true. & ! aggregate + ) & + ] + +end module SnfDislInputModule diff --git a/src/Model/StreamNetworkFlow/snf1flw1.f90 b/src/Model/StreamNetworkFlow/snf1flw1.f90 new file mode 100644 index 00000000000..db937d15dbe --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1flw1.f90 @@ -0,0 +1,357 @@ +!> @brief This module contains the FLW package methods +!! +!! This module can be used to represent inflow to streams. It is based +!! on the GWF WEL package. +!! +!< +module SnfFlwModule + ! -- modules + use KindModule, only: DP, I4B + use ConstantsModule, only: DZERO, DEM1, DONE, LENFTYPE, DNODATA, MAXCHARLEN + use SimVariablesModule, only: errmsg + use SimModule, only: store_error + use MemoryHelperModule, only: create_mem_path + use BndModule, only: BndType + use ObsModule, only: DefaultObsIdProcessor + use SmoothingModule, only: sQSaturation, sQSaturationDerivative + use ObserveModule, only: ObserveType + use TimeSeriesLinkModule, only: TimeSeriesLinkType, & + GetTimeSeriesLinkFromList + use BlockParserModule, only: BlockParserType + use InputOutputModule, only: GetUnit, openfile + use MatrixModule + ! + implicit none + ! + private + public :: flw_create + ! + character(len=LENFTYPE) :: ftype = 'FLW' !< package ftype + character(len=16) :: text = ' FLW' !< package flow text string + ! + type, extends(BndType) :: SnfFlwType + contains + procedure :: bnd_options => flw_options + procedure :: bnd_cf => flw_cf + procedure :: bnd_fc => flw_fc + procedure :: define_listlabel + ! -- methods for observations + procedure, public :: bnd_obs_supported => flw_obs_supported + procedure, public :: bnd_df_obs => flw_df_obs + procedure, public :: bnd_bd_obs => flw_bd_obs + ! -- methods for time series + procedure, public :: bnd_rp_ts => flw_rp_ts + end type SnfFlwType + +contains + + !> @ brief Create a new package object + !! + !! Create a new FLW Package object + !! + !< + subroutine flw_create(packobj, id, ibcnum, inunit, iout, namemodel, pakname) + ! -- dummy variables + class(BndType), pointer :: packobj !< pointer to default package type + integer(I4B), intent(in) :: id !< package id + integer(I4B), intent(in) :: ibcnum !< boundary condition number + integer(I4B), intent(in) :: inunit !< unit number of FLW package input file + integer(I4B), intent(in) :: iout !< unit number of model listing file + character(len=*), intent(in) :: namemodel !< model name + character(len=*), intent(in) :: pakname !< package name + ! -- local variables + type(SnfFlwType), pointer :: flwobj + ! + ! -- allocate the object and assign values to object variables + allocate (flwobj) + packobj => flwobj + ! + ! -- create name and memory path + call packobj%set_names(ibcnum, namemodel, pakname, ftype) + packobj%text = text + ! + ! -- allocate scalars + call flwobj%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 + ! + ! -- return + return + end subroutine flw_create + + !> @ brief Read additional options for package + !! + !! Read additional options for FLW package. + !! + !< + subroutine flw_options(this, option, found) + ! -- modules + use InputOutputModule, only: urword + ! -- dummy variables + class(SnfFlwType), intent(inout) :: this !< SnfFlwType object + character(len=*), intent(inout) :: option !< option keyword string + logical, intent(inout) :: found !< boolean indicating if option found + ! -- local variables + ! -- formats + ! + select case (option) + case ('MOVER') + this%imover = 1 + write (this%iout, '(4x,A)') 'MOVER OPTION ENABLED' + found = .true. + case default + ! + ! -- No options found + found = .false. + end select + ! + ! -- return + return + end subroutine flw_options + + !> @ brief Formulate the package hcof and rhs terms. + !! + !! Formulate the hcof and rhs terms for the FLW package that will be + !! added to the coefficient matrix and right-hand side vector. + !! + !< + subroutine flw_cf(this, reset_mover) + ! -- dummy variables + class(SnfFlwType) :: this !< SnfFlwType object + logical, intent(in), optional :: reset_mover !< boolean for resetting mover + ! -- local variables + integer(I4B) :: i, node + real(DP) :: q + logical :: lrm + ! + ! -- Return if no inflows + if (this%nbound == 0) return + ! + ! -- pakmvrobj cf + lrm = .true. + if (present(reset_mover)) lrm = reset_mover + if (this%imover == 1 .and. lrm) then + call this%pakmvrobj%cf() + end if + ! + ! -- Calculate hcof and rhs for each flw entry + do i = 1, this%nbound + node = this%nodelist(i) + this%hcof(i) = DZERO + if (this%ibound(node) <= 0) then + this%rhs(i) = DZERO + cycle + end if + q = this%bound(1, i) + this%rhs(i) = -q + end do + ! + return + end subroutine flw_cf + + !> @ brief Copy hcof and rhs terms into solution. + !! + !! Add the hcof and rhs terms for the FLW package to the + !! coefficient matrix and right-hand side vector. + !! + !< + subroutine flw_fc(this, rhs, ia, idxglo, matrix_sln) + ! -- dummy variables + class(SnfFlwType) :: this !< SnfFlwType object + real(DP), dimension(:), intent(inout) :: rhs !< right-hand side vector for model + integer(I4B), dimension(:), intent(in) :: ia !< solution CRS row pointers + integer(I4B), dimension(:), intent(in) :: idxglo !< mapping vector for model (local) to solution (global) + class(MatrixBaseType), pointer :: matrix_sln !< solution coefficient matrix + ! -- local variables + integer(I4B) :: i + integer(I4B) :: n + integer(I4B) :: ipos + ! + ! -- pakmvrobj fc + if (this%imover == 1) then + call this%pakmvrobj%fc() + end if + ! + ! -- Copy package rhs and hcof into solution rhs and amat + do i = 1, this%nbound + n = this%nodelist(i) + rhs(n) = rhs(n) + this%rhs(i) + ipos = ia(n) + call matrix_sln%add_value_pos(idxglo(ipos), this%hcof(i)) + ! + ! -- If mover is active and this flw item is discharging, + ! store available water (as positive value). + if (this%imover == 1 .and. this%rhs(i) > DZERO) then + call this%pakmvrobj%accumulate_qformvr(i, this%rhs(i)) + end if + end do + ! + ! -- return + return + end subroutine flw_fc + + !> @ brief Define the list label for the package + !! + !! Method defined the list label for the FLW package. The list label is + !! the heading that is written to iout when PRINT_INPUT option is used. + !! + !< + subroutine define_listlabel(this) + ! -- dummy variables + class(SnfFlwType), intent(inout) :: this !< SnfFlwType object + ! + ! -- create the header list label + this%listlabel = trim(this%filtyp)//' NO.' + if (this%dis%ndim == 3) then + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'ROW' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'COL' + elseif (this%dis%ndim == 2) then + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'LAYER' + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'CELL2D' + else + write (this%listlabel, '(a, a7)') trim(this%listlabel), 'NODE' + end if + write (this%listlabel, '(a, a16)') trim(this%listlabel), 'FLOW RATE' + if (this%inamedbound == 1) then + write (this%listlabel, '(a, a16)') trim(this%listlabel), 'BOUNDARY NAME' + end if + ! + ! -- return + return + end subroutine define_listlabel + + ! -- Procedures related to observations + + !> @brief Determine if observations are supported. + !! + !! Function to determine if observations are supported by the FLW package. + !! Observations are supported by the FLW package. + !! + !! @return flw_obs_supported boolean indicating if observations are supported + !! + !< + logical function flw_obs_supported(this) + ! -- dummy variables + class(SnfFlwType) :: this !< SnfFlwType object + ! + ! -- set boolean + flw_obs_supported = .true. + ! + ! -- return + return + end function flw_obs_supported + + !> @brief Define the observation types available in the package + !! + !! Method to define the observation types available in the FLW package. + !! + !< + subroutine flw_df_obs(this) + ! -- dummy variables + class(SnfFlwType) :: this !< SnfFlwType object + ! -- local variables + integer(I4B) :: indx + ! + ! -- initialize observations + call this%obs%StoreObsType('flw', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor + ! + ! -- Store obs type and assign procedure pointer + ! for to-mvr observation type. + call this%obs%StoreObsType('to-mvr', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => DefaultObsIdProcessor + ! + ! -- return + return + end subroutine flw_df_obs + + !> @brief Save observations for the package + !! + !! Method to save simulated values for the FLW package. + !! + !< + subroutine flw_bd_obs(this) + ! -- dummy variables + class(SnfFlwType) :: this !< SnfFlwType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: n + integer(I4B) :: jj + real(DP) :: v + type(ObserveType), pointer :: obsrv => null() + ! + ! -- clear the observations + call this%obs%obs_bd_clear() + ! + ! -- Save simulated values for all of package's observations. + do i = 1, this%obs%npakobs + obsrv => this%obs%pakobs(i)%obsrv + if (obsrv%BndFound) then + do n = 1, obsrv%indxbnds_count + v = DNODATA + jj = obsrv%indxbnds(n) + select case (obsrv%ObsTypeId) + case ('TO-MVR') + if (this%imover == 1) then + v = this%pakmvrobj%get_qtomvr(jj) + if (v > DZERO) then + v = -v + end if + end if + case ('FLW') + v = this%simvals(jj) + case default + errmsg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId) + call store_error(errmsg) + end select + call this%obs%SaveOneSimval(obsrv, v) + end do + else + call this%obs%SaveOneSimval(obsrv, DNODATA) + end if + end do + ! + ! -- return + return + end subroutine flw_bd_obs + + ! -- Procedure related to time series + + !> @brief Assign time series links for the package + !! + !! Assign the time series links for the FLW package. Only + !! the Q variable can be defined with time series. + !! + !< + subroutine flw_rp_ts(this) + ! -- dummy variables + class(SnfFlwType), intent(inout) :: this !< SnfFlwType object + ! -- local variables + integer(I4B) :: i, nlinks + type(TimeSeriesLinkType), pointer :: tslink => null() + ! + ! -- set up the time series links + nlinks = this%TsManager%boundtslinks%Count() + do i = 1, nlinks + tslink => GetTimeSeriesLinkFromList(this%TsManager%boundtslinks, i) + if (associated(tslink)) then + if (tslink%JCol == 1) then + tslink%Text = 'Q' + end if + end if + end do + ! + ! -- return + return + end subroutine flw_rp_ts + +end module SnfFlwModule diff --git a/src/Model/StreamNetworkFlow/snf1mmr1.f90 b/src/Model/StreamNetworkFlow/snf1mmr1.f90 new file mode 100644 index 00000000000..47420aae58e --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1mmr1.f90 @@ -0,0 +1,978 @@ +!> @brief Stream Network Flow (SNF) Muskingum-Manning Routing (MMR) Module +!! +!! This module solves one-dimension flow routing using a Muskingum-Manning +!! approach. The approach uses the following input parameters. +!! +!< +module SnfMmrModule + + use KindModule, only: DP, I4B, LGP + use ConstantsModule, only: LENMEMPATH, LENVARNAME, LINELENGTH, & + DZERO, DHALF, DONE, DTWO, DTHREE, & + LENBUDTXT + use MemoryHelperModule, only: create_mem_path + use MemoryManagerModule, only: mem_allocate + use SimVariablesModule, only: errmsg, warnmsg + use SimModule, only: count_errors, store_error, store_error_unit + use NumericalPackageModule, only: NumericalPackageType + use BaseDisModule, only: DisBaseType + use SnfDislModule, only: SnfDislType + use ObsModule, only: ObsType, obs_cr + use ObsModule, only: DefaultObsIdProcessor + use ObserveModule, only: ObserveType + + implicit none + private + public :: SnfMmrType, mmr_cr + + + character(len=LENBUDTXT), dimension(2) :: budtxt = & !< text labels for budget terms + &[' STORAGE', ' EXT-OUTFLOW'] + + type, extends(NumericalPackageType) :: SnfMmrType + + ! -- user-provided input + integer(I4B), dimension(:), pointer, contiguous :: iseg_order => null() !< routing calculation order + real(DP), dimension(:), pointer, contiguous :: qoutflow0 => null() !< initial outflow for each reach + real(DP), dimension(:), pointer, contiguous :: k_coef => null() !< manning K coefficient + real(DP), dimension(:), pointer, contiguous :: x_coef => null() !< routing weighting factor + + ! -- input arguments to calc_muskingum_mann routine + real(DP), dimension(:), pointer, contiguous :: inflow_old => null() !< inflow to each segment for last time step + real(DP), dimension(:), pointer, contiguous :: inflow_new => null() !< inflow to each segment for current time step + real(DP), dimension(:), pointer, contiguous :: outflow_old => null() !< outflow from each segment for last time step + real(DP), dimension(:), pointer, contiguous :: outflow_new => null() !< outflow from each segment for current time step + real(DP), dimension(:), pointer, contiguous :: c0 => null() !< Muskingum c0 variable + real(DP), dimension(:), pointer, contiguous :: c1 => null() !< Muskingum c1 variable + real(DP), dimension(:), pointer, contiguous :: c2 => null() !< Muskingum c2 variable + + ! -- budget vectors + real(DP), dimension(:), pointer, contiguous :: qextoutflow => null() !< flows leaving model (for tosegment = 0) + real(DP), dimension(:), pointer, contiguous :: qsto => null() !< storage rates + + ! -- pointer to concrete disl subclass of DisBaseType + type(SnfDislType), pointer :: disl + + ! -- observation data + integer(I4B), pointer :: inobspkg => null() !< unit number for obs package + type(ObsType), pointer :: obs => null() !< observation package + + contains + + procedure :: allocate_scalars + procedure :: allocate_arrays + procedure :: mmr_load + procedure :: source_options + procedure :: log_options + procedure :: source_griddata + procedure :: log_griddata + procedure :: mmr_ar + procedure :: mmr_rp + procedure :: mmr_ad + procedure :: mmr_init_data + procedure :: mmr_solve + procedure :: mmr_cq + procedure :: mmr_bd + procedure :: mmr_save_model_flows + procedure :: mmr_print_model_flows + procedure :: mmr_da + procedure :: mmr_df_obs + procedure :: mmr_rp_obs + procedure :: mmr_bd_obs + + end type SnfMmrType + + contains + + !> @brief create package + !< + subroutine mmr_cr(mmrobj, name_model, inunit, dis, iout) + ! -- modules + use IdmMf6FileLoaderModule, only: input_load + use ConstantsModule, only: LENPACKAGETYPE + ! -- dummy + type(SnfMmrType), pointer :: mmrobj + character(len=*), intent(in) :: name_model + integer(I4B), intent(in) :: inunit + class(DisBaseType), pointer, intent(inout) :: dis !< the pointer to the discretization + integer(I4B), intent(in) :: iout + ! -- formats + character(len=*), parameter :: fmtheader = & + "(1x, /1x, 'MMR -- MUSKINGUM MANNINGS ROUTING PACKAGE, VERSION 1, 1/23/2023', & + &' INPUT READ FROM UNIT ', i0, /)" + ! + ! -- Create the object + allocate (mmrobj) + call mmrobj%set_names(1, name_model, 'MMR', 'MMR') + call mmrobj%allocate_scalars() + mmrobj%inunit = inunit + mmrobj%iout = iout + mmrobj%dis => dis + + ! -- store pointer to disl + ! Not normally good practice, but since SNF only works with DISL + ! may be okay + select type (dis) + type is (SnfDislType) + mmrobj%disl => dis + end select + + ! -- create obs package + call obs_cr(mmrobj%obs, mmrobj%inobspkg) + + ! + ! -- if reading from file + if (inunit > 0) then + ! + ! -- Identify package + if (iout > 0) then + write (iout, fmtheader) inunit + end if + ! + ! -- allocate arrays + call mmrobj%allocate_arrays() + ! + ! -- load mmr + call mmrobj%mmr_load() + end if + ! + ! -- Return + return + end subroutine mmr_cr + + !> @ brief Allocate scalars + !! + !! Allocate and initialize scalars for the package. The base model + !! allocate scalars method is also called. + !! + !< + subroutine allocate_scalars(this) + ! -- modules + ! -- dummy + class(SnfMmrtype) :: this + ! + ! -- allocate scalars in NumericalPackageType + call this%NumericalPackageType%allocate_scalars() + ! + ! -- Allocate scalars + call mem_allocate(this%inobspkg, 'INOBSPKG', this%memoryPath) + + this%inobspkg = 0 + + return + end subroutine allocate_scalars + + !> @brief allocate memory for arrays + !< + subroutine allocate_arrays(this) + ! -- dummy + class(SnfMmrType) :: this + ! -- locals + integer(I4B) :: n + ! + ! -- user-provided input + call mem_allocate(this%iseg_order, this%dis%nodes, 'ISEG_ORDER', this%memoryPath) + call mem_allocate(this%qoutflow0, this%dis%nodes, 'QOUTFLOW0', this%memoryPath) + call mem_allocate(this%k_coef, this%dis%nodes, 'K_COEF', this%memoryPath) + call mem_allocate(this%x_coef, this%dis%nodes, 'X_COEF', this%memoryPath) + + ! -- input arguments to calc_muskingum_mann routine + call mem_allocate(this%inflow_old, this%dis%nodes, 'INFLOW_OLD', this%memoryPath) + call mem_allocate(this%inflow_new, this%dis%nodes, 'INFLOW_NEW', this%memoryPath) + call mem_allocate(this%outflow_old, this%dis%nodes, 'OUTFLOW_OLD', this%memoryPath) + call mem_allocate(this%outflow_new, this%dis%nodes, 'OUTFLOW_NEW', this%memoryPath) + call mem_allocate(this%c0, this%dis%nodes, 'C0', this%memoryPath) + call mem_allocate(this%c1, this%dis%nodes, 'C1', this%memoryPath) + call mem_allocate(this%c2, this%dis%nodes, 'C2', this%memoryPath) + + ! -- budgeting variables + call mem_allocate(this%qextoutflow, this%dis%nodes, 'QEXTOUTFLOW', this%memoryPath) + call mem_allocate(this%qsto, this%dis%nodes, 'QSTO', this%memoryPath) + + do n = 1, this%dis%nodes + + this%iseg_order(n) = 0 + this%qoutflow0(n) = DZERO + this%k_coef(n) = DZERO + this%x_coef(n) = DZERO + + this%inflow_old(n) = DZERO + this%inflow_new(n) = DZERO + this%outflow_old(n) = DZERO + this%outflow_new(n) = DZERO + this%c0(n) = DZERO + this%c1(n) = DZERO + this%c2(n) = DZERO + + this%qextoutflow(n) = DZERO + this%qsto(n) = DZERO + + end do + + ! -- Return + return + end subroutine allocate_arrays + + !> @brief load data from IDM to package + !< + subroutine mmr_load(this) + ! -- dummy + class(SnfMmrType) :: this + ! -- locals + ! + ! -- source input data + call this%source_options() + call this%source_griddata() + ! + ! -- Return + return + end subroutine mmr_load + + !> @brief Copy options from IDM into package + !< + subroutine source_options(this) + ! -- modules + use KindModule, only: LGP + use InputOutputModule, only: getunit, openfile + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + use SnfMmrInputModule, only: SnfMmrParamFoundType + ! -- dummy + class(SnfMmrType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + type(SnfMmrParamFoundType) :: found + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'MMR', idm_context) + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%iprflow, 'IPRFLOW', idmMemoryPath, found%iprflow) + call mem_set_value(this%ipakcb, 'IPAKCB', idmMemoryPath, found%ipakcb) + call mem_set_value(this%obs%inputFilename, 'OBS6_FILENAME', idmMemoryPath, & + found%obs6_filename) + ! + ! -- save flows option active + if (found%ipakcb) this%ipakcb = -1 + + if (found%obs6_filename) then + this%obs%active = .true. + this%inobspkg = GetUnit() + this%obs%inUnitObs = this%inobspkg + call openfile(this%inobspkg, this%iout, this%obs%inputFilename, 'OBS') + call this%obs%obs_df(this%iout, this%packName, this%filtyp, this%dis) + call this%mmr_df_obs() + end if + + ! + ! -- log values to list file + if (this%iout > 0) then + call this%log_options(found) + end if + ! + ! -- Return + return + end subroutine source_options + + !> @brief Write user options to list file + !< + subroutine log_options(this, found) + use SnfMmrInputModule, only: SnfMmrParamFoundType + class(SnfMmrType) :: this + type(SnfMmrParamFoundType), intent(in) :: found + + write (this%iout, '(1x,a)') 'Setting MMR Options' + + if (found%iprflow) then + write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed & + &to listing file whenever ICBCFL is not zero.' + end if + + if (found%ipakcb) then + write (this%iout, '(4x,a)') 'Cell-by-cell flow information will be printed & + &to listing file whenever ICBCFL is not zero.' + end if + + + write (this%iout, '(1x,a,/)') 'End Setting MMR Options' + + end subroutine log_options + + !> @brief copy griddata from IDM to package + !< + subroutine source_griddata(this) + ! -- modules + use SimModule, only: count_errors, store_error + use MemoryHelperModule, only: create_mem_path + use MemoryManagerModule, only: mem_reallocate + use MemoryManagerExtModule, only: mem_set_value + use SimVariablesModule, only: idm_context + use SnfMmrInputModule, only: SnfMmrParamFoundType + ! -- dummy + class(SnfMmrType) :: this + ! -- locals + character(len=LENMEMPATH) :: idmMemoryPath + type(SnfMmrParamFoundType) :: found + integer(I4B), dimension(:), pointer, contiguous :: map + ! -- formats + ! + ! -- set memory path + idmMemoryPath = create_mem_path(this%name_model, 'MMR', idm_context) + ! + ! -- set map to convert user input data into reduced data + map => null() + if (this%dis%nodes < this%dis%nodesuser) map => this%dis%nodeuser + ! + ! -- update defaults with idm sourced values + call mem_set_value(this%iseg_order, 'ISEG_ORDER', idmMemoryPath, map, found%iseg_order) + call mem_set_value(this%qoutflow0, 'QOUTFLOW0', idmMemoryPath, map, found%qoutflow0) + call mem_set_value(this%k_coef, 'K_COEF', idmMemoryPath, map, found%k_coef) + call mem_set_value(this%x_coef, 'X_COEF', idmMemoryPath, map, found%x_coef) + ! + ! -- ensure ISEG_ORDER was found + if (.not. found%iseg_order) then + write (errmsg, '(a)') 'Error in GRIDDATA block: ISEG_ORDER not found.' + call store_error(errmsg) + end if + ! + ! -- ensure QOUTFLOW0 was found + if (.not. found%qoutflow0) then + write (errmsg, '(a)') 'Error in GRIDDATA block: QOUTFLOW0 not found.' + call store_error(errmsg) + end if + ! + ! -- ensure MANN_N was found + if (.not. found%k_coef) then + write (errmsg, '(a)') 'Error in GRIDDATA block: K_COEF not found.' + call store_error(errmsg) + end if + ! + ! -- ensure X_COEF was found + if (.not. found%x_coef) then + write (errmsg, '(a)') 'Error in GRIDDATA block: X_COEF not found.' + call store_error(errmsg) + end if + ! + ! -- log griddata + if (this%iout > 0) then + call this%log_griddata(found) + end if + ! + ! -- Return + return + end subroutine source_griddata + + !> @brief log griddata to list file + !< + subroutine log_griddata(this, found) + use SnfMmrInputModule, only: SnfMmrParamFoundType + class(SnfMmrType) :: this + type(SnfMmrParamFoundType), intent(in) :: found + + write (this%iout, '(1x,a)') 'Setting MMR Griddata' + + if (found%iseg_order) then + write (this%iout, '(4x,a)') 'ISEG_ORDER set from input file' + end if + + if (found%qoutflow0) then + write (this%iout, '(4x,a)') 'QOUTFLOW0 set from input file' + end if + + if (found%k_coef) then + write (this%iout, '(4x,a)') 'K_COEF set from input file' + end if + + if (found%x_coef) then + write (this%iout, '(4x,a)') 'X_COEF set from input file' + end if + + write (this%iout, '(1x,a,/)') 'End Setting MMR Griddata' + + end subroutine log_griddata + + !> @brief allocate memory + !< + subroutine mmr_ar(this) + ! -- modules + ! -- dummy + class(SnfMmrType) :: this !< this instance + ! + + ! - observation data + call this%obs%obs_ar() + + ! -- initialize routing variables + call this%mmr_init_data() + + return + end subroutine mmr_ar + + !> @brief allocate memory + !< + subroutine mmr_rp(this) + ! -- modules + ! -- dummy + class(SnfMmrType) :: this !< this instance + ! + ! -- read observations + call this%mmr_rp_obs() + return + end subroutine mmr_rp + + subroutine mmr_ad(this, irestore) + ! + class(SnfMmrType) :: this + integer(I4B), intent(in) :: irestore + integer(I4B) :: n + ! + ! Advance forward or backward depending on irestore + if (irestore == 0) then + do n = 1, this%disl%nodes + this%inflow_old(n) = this%inflow_new(n) + this%outflow_old(n) = this%outflow_new(n) + end do + else + do n = 1, this%disl%nodes + this%inflow_new(n) = this%inflow_old(n) + this%outflow_new(n) = this%outflow_old(n) + end do + end if + + ! -- Push simulated values to preceding time/subtime step + call this%obs%obs_ad() + + ! + ! -- Return + return + end subroutine mmr_ad + + !> @brief solve + !< + subroutine mmr_solve(this, rhs) + ! -- modules + use TdisModule, only: delt + ! -- dummy + class(SnfMmrType) :: this !< this instance + real(DP), dimension(:), intent(in) :: rhs !< right-hand-side vector of boundary package inflows + ! -- local + + call calc_muskingum_coefficients(delt, this%k_coef, this%x_coef, this%c0, & + this%c1, this%c2) + + call calc_muskingum(this%disl%tosegment, this%iseg_order, this%inflow_old, & + this%outflow_old, this%inflow_new, this%outflow_new, & + this%c0, this%c1, this%c2, -rhs) + + ! -- return + return + end subroutine mmr_solve + + subroutine mmr_cq(this, flowja) + ! -- dummy + class(SnfMmrType) :: this + real(DP), intent(inout), dimension(:) :: flowja + ! -- local + integer(I4B) :: n, ipos, m + real(DP) :: qnm, q + ! + ! -- Transfer seg_outflow into flowja + do n = 1, this%dis%nodes + m = this%disl%tosegment(n) + ! TODO: may be a faster way than lookup + if (m > 0) then + do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 + if (this%dis%con%ja(ipos) == m) exit + end do + qnm = this%outflow_new(n) + flowja(ipos) = -qnm + flowja(this%dis%con%isym(ipos)) = qnm + end if + end do + + ! Transfer any flows leaving tosegment 0 into qextoutflow + do n = 1, this%dis%nodes + q = DZERO + if (this%disl%tosegment(n) == 0) then + q = -this%outflow_new(n) + end if + this%qextoutflow(n) = q + ! + ! -- add to diagonal of flowja + ipos = this%dis%con%ia(n) + flowja(ipos) = flowja(ipos) + q + end do + + ! Transfer storage terms into qsto + do n = 1, this%dis%nodes + q = this%outflow_new(n) - this%inflow_new(n) + this%qsto(n) = q + ! + ! -- add to diagonal + ipos = this%dis%con%ia(n) + flowja(ipos) = flowja(ipos) + q + end do + + ! + ! -- Return + return + end subroutine mmr_cq + + !> @ brief Model budget calculation for package + !! + !! Budget calculation for the MMR package components. Components include + !! external outflow + !! + !< + subroutine mmr_bd(this, isuppress_output, model_budget) + ! -- modules + use TdisModule, only: delt + use BudgetModule, only: BudgetType, rate_accumulator + ! -- dummy variables + class(SnfMmrType) :: this !< SnfMmrType object + integer(I4B), intent(in) :: isuppress_output !< flag to suppress model output + type(BudgetType), intent(inout) :: model_budget !< model budget object + ! -- local variables + real(DP) :: rin + real(DP) :: rout + ! + ! -- Add external outflow rates to model budget + call rate_accumulator(this%qextoutflow, rin, rout) + call model_budget%addentry(rin, rout, delt, ' MMR', & + isuppress_output, ' EXT-OUTFLOW') + ! + ! -- Add storage rates to model budget + call rate_accumulator(this%qsto, rin, rout) + call model_budget%addentry(rin, rout, delt, ' MMR', & + isuppress_output, ' STORAGE') + ! + ! -- return + return + end subroutine mmr_bd + + !> @ brief save flows for package + !< + subroutine mmr_save_model_flows(this, flowja, icbcfl, icbcun) + ! -- dummy + class(SnfMmrType) :: this + real(DP), dimension(:), intent(in) :: flowja + integer(I4B), intent(in) :: icbcfl + integer(I4B), intent(in) :: icbcun + ! -- local + real(DP) :: dinact + integer(I4B) :: ibinun + integer(I4B) :: iprint, nvaluesp, nwidthp + character(len=1) :: cdatafmp = ' ', editdesc = ' ' + ! -- formats + ! + ! -- Set unit number for binary output + iprint = 0 + dinact = DZERO + if (this%ipakcb < 0) then + ibinun = icbcun + elseif (this%ipakcb == 0) then + ibinun = 0 + else + ibinun = this%ipakcb + end if + if (icbcfl == 0) ibinun = 0 + ! + ! -- Write the face flows if requested + if (ibinun /= 0) then + ! + ! -- flowja + call this%dis%record_connection_array(flowja, ibinun, this%iout) + ! + ! + ! -- storage + call this%dis%record_array(this%qsto, this%iout, iprint, -ibinun, & + budtxt(1), cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + + ! -- external outflow + ! TODO: should this be written as a list instead? + call this%dis%record_array(this%qextoutflow, this%iout, iprint, -ibinun, & + budtxt(2), cdatafmp, nvaluesp, & + nwidthp, editdesc, dinact) + + end if + ! + ! -- Return + return + end subroutine mmr_save_model_flows + + !> @ brief print flows for package + !< + subroutine mmr_print_model_flows(this, ibudfl, flowja) + ! -- modules + use TdisModule, only: kper, kstp + use ConstantsModule, only: LENBIGLINE + ! -- dummy + class(SnfMmrType) :: this + integer(I4B), intent(in) :: ibudfl + real(DP), intent(inout), dimension(:) :: flowja + ! -- local + character(len=LENBIGLINE) :: line + character(len=30) :: tempstr + integer(I4B) :: n, ipos, m + real(DP) :: qnm + ! -- formats + character(len=*), parameter :: fmtiprflow = & + &"(/,4x,'CALCULATED INTERCELL FLOW FOR PERIOD ', i0, ' STEP ', i0)" +! ------------------------------------------------------------------------------ + ! + ! -- Write flowja to list file if requested + if (ibudfl /= 0 .and. this%iprflow > 0) then + write (this%iout, fmtiprflow) kper, kstp + do n = 1, this%dis%nodes + line = '' + call this%dis%noder_to_string(n, tempstr) + line = trim(tempstr)//':' + do ipos = this%dis%con%ia(n) + 1, this%dis%con%ia(n + 1) - 1 + m = this%dis%con%ja(ipos) + call this%dis%noder_to_string(m, tempstr) + line = trim(line)//' '//trim(tempstr) + qnm = flowja(ipos) + write (tempstr, '(1pg15.6)') qnm + line = trim(line)//' '//trim(adjustl(tempstr)) + end do + write (this%iout, '(a)') trim(line) + end do + end if + ! + ! -- Return + return + end subroutine mmr_print_model_flows + + !> @brief deallocate memory + !< + subroutine mmr_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + use MemoryManagerExtModule, only: memorylist_remove + use SimVariablesModule, only: idm_context + ! -- dummy + class(SnfMmrType) :: this + ! + ! -- Deallocate input memory + call memorylist_remove(this%name_model, 'MMR', idm_context) + ! + ! -- Scalars + ! + ! -- Deallocate arrays + call mem_deallocate(this%iseg_order) + call mem_deallocate(this%qoutflow0) + call mem_deallocate(this%k_coef) + call mem_deallocate(this%x_coef) + + ! -- input arguments to calc_muskingum_mann routine + call mem_deallocate(this%inflow_old) + call mem_deallocate(this%inflow_new) + call mem_deallocate(this%outflow_old) + call mem_deallocate(this%outflow_new) + call mem_deallocate(this%c0) + call mem_deallocate(this%c1) + call mem_deallocate(this%c2) + + ! -- budget variables + call mem_deallocate(this%qextoutflow) + call mem_deallocate(this%qsto) + + ! -- obs package + call mem_deallocate(this%inobspkg) + call this%obs%obs_da() + deallocate (this%obs) + nullify (this%obs) + + ! -- deallocate parent + call this%NumericalPackageType%da() + ! + ! -- Return + return + end subroutine mmr_da + + !> @brief initialize mmr data + !! + !! This routine is only called once at the beginning of + !! the simulation from mmr_ar() + !! + !< + subroutine mmr_init_data(this) + ! -- modules + ! -- dummy + class(SnfMmrType) :: this !< this instance + ! -- local + integer(I4B) :: n + integer(I4B) :: ito_seg + integer(I4B), dimension(:), allocatable :: nupreaches + + ! -- Set inflow_new and outflow_new from qoutflow0 + ! These will be copied into inflow_old and outflow_old + ! as part of the first advance + do n = 1, this%dis%nodes + this%outflow_new(n) = this%qoutflow0(n) + ito_seg = this%disl%tosegment(n) + if (ito_seg > 0) then + this%inflow_new(ito_seg) = & + this%inflow_new(ito_seg) + this%qoutflow0(n) + end if + end do + + ! -- Calculate the number of upstream reaches connected to each reach, + ! which is then used in the subsequent section to initialize the + ! flow conditions + allocate(nupreaches(this%dis%nodes)) + do n = 1, this%dis%nodes + nupreaches(n) = 0 + end do + do n = 1, this%dis%nodes + ito_seg = this%disl%tosegment(n) + if (ito_seg > 0) then + nupreaches(ito_seg) = nupreaches(ito_seg) + 1 + end if + end do + + ! -- For any cells that do not have inflow via upstream reaches, + ! the inflow_new term should be set to the initial outflow + ! for the reach, qoutflow0 + do n = 1, this%dis%nodes + if (nupreaches(n) == 0) then + this%inflow_new(n) = this%qoutflow0(n) + end if + end do + + ! -- return + return + end subroutine mmr_init_data + + !> @brief calculate muskingum c0, c1, c2 coefficients + !< + subroutine calc_muskingum_coefficients(delt, k_coef, x_coef, c0, c1, c2) + ! -- modules + ! -- dummy + real(DP), intent(in) :: delt !< Length of time step + real(DP), dimension(:), intent(in) :: k_coef !< Muskingum K coefficient + real(DP), dimension(:), intent(in) :: x_coef !< Muskingum K coefficient + real(DP), dimension(:), intent(inout) :: c0 !< Muskingum K coefficient + real(DP), dimension(:), intent(inout) :: c1 !< Muskingum K coefficient + real(DP), dimension(:), intent(inout) :: c2 !< Muskingum K coefficient + ! -- local + integer(I4B) :: n + real(DP) :: dtoverk + real(DP) :: twox + real(DP) :: two_oneminusx + real(DP) :: denom + + ! -- Calculate muskingum C coefficients, using Ponce equations + do n = 1, size(k_coef) + dtoverk = delt / k_coef(n) + twox = DTWO * x_coef(n) + two_oneminusx = DTWO * (DONE - x_coef(n)) + denom = two_oneminusx + dtoverk + c0(n) = (dtoverk - twox) / denom + c1(n) = (dtoverk + twox) / denom + c2(n) = (two_oneminusx - dtoverk) / denom + end do + + ! -- return + return + end subroutine calc_muskingum_coefficients + + subroutine calc_muskingum(itosegment, iseg_order, inflow_old, outflow_old, & + inflow_new, outflow_new, c0, c1, c2, qsource) + ! -- dummy + integer(I4B), dimension(:), intent(in) :: itosegment + integer(I4B), dimension(:), intent(in) :: iseg_order + real(DP), dimension(:), intent(in) :: inflow_old + real(DP), dimension(:), intent(in) :: outflow_old + real(DP), dimension(:), intent(inout) :: inflow_new + real(DP), dimension(:), intent(inout) :: outflow_new + real(DP), dimension(:), intent(in) :: c0 + real(DP), dimension(:), intent(in) :: c1 + real(DP), dimension(:), intent(in) :: c2 + real(DP), dimension(:), intent(in) :: qsource + ! -- dummy + integer(I4B) :: n + integer(I4B) :: j + integer(I4B) :: i + real(DP) :: qoutflow + real(DP) :: qinflow + + ! -- Initialize inflow_new with any sources, such as lateral inflow + do n = 1, size(itosegment) + j = iseg_order(n) + inflow_new(j) = qsource(j) + end do + + ! -- Use Muskingum method to calculate outflows and accumulate + ! outflows in downstream reaches + do n = 1, size(itosegment) + j = iseg_order(n) + qinflow = inflow_new(j) !+ qsource(j) + qoutflow = c0(j) * qinflow + c1(j) * inflow_old(j) + c2(j) * outflow_old(j) + outflow_new(j) = qoutflow + i = itosegment(j) + if (i > 0) then + inflow_new(i) = inflow_new(i) + qoutflow + end if + end do + + end subroutine calc_muskingum + + !> @brief Define the observation types available in the package + !! + !! Method to define the observation types available in the package. + !! + !< + subroutine mmr_df_obs(this) + ! -- dummy variables + class(SnfMmrType) :: this !< SnfMmrType object + ! -- local variables + integer(I4B) :: indx + ! + ! -- Store obs type and assign procedure pointer + ! for ext-outflow observation type. + call this%obs%StoreObsType('ext-outflow', .true., indx) + this%obs%obsData(indx)%ProcessIdPtr => mmrobsidprocessor + ! + ! -- return + return + end subroutine mmr_df_obs + + + subroutine mmrobsidprocessor(obsrv, dis, inunitobs, iout) + ! -- dummy + type(ObserveType), intent(inout) :: obsrv + class(DisBaseType), intent(in) :: dis + integer(I4B), intent(in) :: inunitobs + integer(I4B), intent(in) :: iout + ! -- local + integer(I4B) :: n + character(len=LINELENGTH) :: strng + ! + ! -- Initialize variables + strng = obsrv%IDstring + read(strng, *) n + ! + if (n > 0) then + obsrv%NodeNumber = n + else + errmsg = 'Error reading data from ID string' + call store_error(errmsg) + call store_error_unit(inunitobs) + end if + ! + return + end subroutine mmrobsidprocessor + + + !> @brief Save observations for the package + !! + !! Method to save simulated values for the package. + !! + !< + subroutine mmr_bd_obs(this) + ! -- dummy variables + class(SnfMmrType) :: this !< SnfMmrType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: j + integer(I4B) :: n + real(DP) :: v + character(len=100) :: msg + type(ObserveType), pointer :: obsrv => null() + ! + ! Write simulated values for all observations + if (this%obs%npakobs > 0) then + call this%obs%obs_bd_clear() + do i = 1, this%obs%npakobs + obsrv => this%obs%pakobs(i)%obsrv + do j = 1, obsrv%indxbnds_count + n = obsrv%indxbnds(j) + v = DZERO + select case (obsrv%ObsTypeId) + case ('EXT-OUTFLOW') + v = this%qextoutflow(n) + case default + msg = 'Unrecognized observation type: '//trim(obsrv%ObsTypeId) + call store_error(msg) + end select + call this%obs%SaveOneSimval(obsrv, v) + end do + end do + ! + ! -- write summary of package error messages + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + end if + ! + ! -- return + return + end subroutine mmr_bd_obs + + !> @brief Read and prepare observations for a package + !! + !! Method to read and prepare observations for a package. + !! + !< + subroutine mmr_rp_obs(this) + ! -- modules + use TdisModule, only: kper + ! -- dummy variables + class(SnfMmrType), intent(inout) :: this !< SnfMmrType object + ! -- local variables + integer(I4B) :: i + integer(I4B) :: j + integer(I4B) :: nn1 + class(ObserveType), pointer :: obsrv => null() + ! -- formats + ! + ! -- process each package observation + ! only done the first stress period since boundaries are fixed + ! for the simulation + if (kper == 1) then + do i = 1, this%obs%npakobs + obsrv => this%obs%pakobs(i)%obsrv + ! + ! -- get node number 1 + nn1 = obsrv%NodeNumber + if (nn1 < 1 .or. nn1 > this%disl%nodes) then + write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') & + trim(adjustl(obsrv%ObsTypeId)), & + 'reach must be greater than 0 and less than or equal to', & + this%disl%nodes, '(specified value is ', nn1, ')' + call store_error(errmsg) + else + if (obsrv%indxbnds_count == 0) then + call obsrv%AddObsIndex(nn1) + else + errmsg = 'Programming error in mmr_rp_obs' + call store_error(errmsg) + end if + end if + ! + ! -- check that node number 1 is valid; call store_error if not + do j = 1, obsrv%indxbnds_count + nn1 = obsrv%indxbnds(j) + if (nn1 < 1 .or. nn1 > this%disl%nodes) then + write (errmsg, '(a,1x,a,1x,i0,1x,a,1x,i0,a)') & + trim(adjustl(obsrv%ObsTypeId)), & + 'reach must be greater than 0 and less than or equal to', & + this%disl%nodes, '(specified value is ', nn1, ')' + call store_error(errmsg) + end if + end do + end do + ! + ! -- evaluate if there are any observation errors + if (count_errors() > 0) then + call this%parser%StoreErrorUnit() + end if + end if + ! + ! -- return + return + end subroutine mmr_rp_obs + +end module SnfMmrModule \ No newline at end of file diff --git a/src/Model/StreamNetworkFlow/snf1mmr1idm.f90 b/src/Model/StreamNetworkFlow/snf1mmr1idm.f90 new file mode 100644 index 00000000000..d1ea475aa74 --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1mmr1idm.f90 @@ -0,0 +1,201 @@ +module SnfMmrInputModule + use InputDefinitionModule, only: InputParamDefinitionType, & + InputBlockDefinitionType + private + public snf_mmr_param_definitions + public snf_mmr_aggregate_definitions + public snf_mmr_block_definitions + public SnfMmrParamFoundType + + type SnfMmrParamFoundType + logical :: ipakcb = .false. + logical :: iprflow = .false. + logical :: obs_filerecord = .false. + logical :: obs6 = .false. + logical :: obs6_filename = .false. + logical :: iseg_order = .false. + logical :: qoutflow0 = .false. + logical :: k_coef = .false. + logical :: x_coef = .false. + end type SnfMmrParamFoundType + + type(InputParamDefinitionType), parameter :: & + snfmmr_ipakcb = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'OPTIONS', & ! block + 'SAVE_FLOWS', & ! tag name + 'IPAKCB', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_iprflow = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'OPTIONS', & ! block + 'PRINT_FLOWS', & ! tag name + 'IPRFLOW', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_obs_filerecord = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'OPTIONS', & ! block + 'OBS_FILERECORD', & ! tag name + 'OBS_FILERECORD', & ! fortran variable + 'RECORD OBS6 FILEIN OBS6_FILENAME', & ! type + '', & ! shape + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_obs6 = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6', & ! tag name + 'OBS6', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_obs6_filename = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'OPTIONS', & ! block + 'OBS6_FILENAME', & ! tag name + 'OBS6_FILENAME', & ! fortran variable + 'STRING', & ! type + '', & ! shape + .true., & ! required + .true., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_iseg_order = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'GRIDDATA', & ! block + 'ISEG_ORDER', & ! tag name + 'ISEG_ORDER', & ! fortran variable + 'INTEGER1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_qoutflow0 = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'GRIDDATA', & ! block + 'QOUTFLOW0', & ! tag name + 'QOUTFLOW0', & ! fortran variable + 'DOUBLE1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_k_coef = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'GRIDDATA', & ! block + 'K_COEF', & ! tag name + 'K_COEF', & ! fortran variable + 'DOUBLE1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snfmmr_x_coef = InputParamDefinitionType & + ( & + 'SNF', & ! component + 'MMR', & ! subcomponent + 'GRIDDATA', & ! block + 'X_COEF', & ! tag name + 'X_COEF', & ! fortran variable + 'DOUBLE1D', & ! type + 'NODES', & ! shape + .true., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false. & ! layered + ) + + type(InputParamDefinitionType), parameter :: & + snf_mmr_param_definitions(*) = & + [ & + snfmmr_ipakcb, & + snfmmr_iprflow, & + snfmmr_obs_filerecord, & + snfmmr_obs6, & + snfmmr_obs6_filename, & + snfmmr_iseg_order, & + snfmmr_qoutflow0, & + snfmmr_k_coef, & + snfmmr_x_coef & + ] + + type(InputParamDefinitionType), parameter :: & + snf_mmr_aggregate_definitions(*) = & + [ & + InputParamDefinitionType :: & + ] + + type(InputBlockDefinitionType), parameter :: & + snf_mmr_block_definitions(*) = & + [ & + InputBlockDefinitionType( & + 'OPTIONS', & ! blockname + .false., & ! required + .false. & ! aggregate + ), & + InputBlockDefinitionType( & + 'GRIDDATA', & ! blockname + .true., & ! required + .false. & ! aggregate + ) & + ] + +end module SnfMmrInputModule diff --git a/src/Model/StreamNetworkFlow/snf1oc1.f90 b/src/Model/StreamNetworkFlow/snf1oc1.f90 new file mode 100644 index 00000000000..26c647c4b1e --- /dev/null +++ b/src/Model/StreamNetworkFlow/snf1oc1.f90 @@ -0,0 +1,101 @@ +module SnfOcModule + + use BaseDisModule, only: DisBaseType + use KindModule, only: DP, I4B + use ConstantsModule, only: LENMODELNAME + use OutputControlModule, only: OutputControlType + use OutputControlDataModule, only: OutputControlDataType, ocd_cr + + implicit none + private + public SnfOcType, oc_cr + + !> @ brief Output control for GWF + !! + !! Concrete implementation of OutputControlType for the + !! GWF Model + !< + type, extends(OutputControlType) :: SnfOcType + contains + procedure :: oc_ar + end type SnfOcType + +contains + + !> @ brief Create SnfOcType + !! + !! Create by allocating a new SnfOcType object and initializing + !! member variables. + !! + !< + subroutine oc_cr(ocobj, name_model, inunit, iout) + ! -- dummy + type(SnfOcType), pointer :: ocobj !< SnfOcType object + character(len=*), intent(in) :: name_model !< name of the model + integer(I4B), intent(in) :: inunit !< unit number for input + integer(I4B), intent(in) :: iout !< unit number for output + ! + ! -- Create the object + allocate (ocobj) + ! + ! -- Allocate scalars + call ocobj%allocate_scalars(name_model) + ! + ! -- Save unit numbers + ocobj%inunit = inunit + ocobj%iout = iout + ! + ! -- Initialize block parser + call ocobj%parser%Initialize(inunit, iout) + ! + ! -- Return + return + end subroutine oc_cr + + !> @ brief Allocate and read SnfOcType + !! + !! Setup head and budget as output control variables. + !! + !< + subroutine oc_ar(this, stage, dis, dnodata) + ! -- dummy + class(SnfOcType) :: this !< SnfOcType object + real(DP), dimension(:), pointer, contiguous, intent(in) :: stage !< model stage + class(DisBaseType), pointer, intent(in) :: dis !< model discretization package + real(DP), intent(in) :: dnodata !< no data value + ! -- local + integer(I4B) :: i, nocdobj, inodata + type(OutputControlDataType), pointer :: ocdobjptr + real(DP), dimension(:), pointer, contiguous :: nullvec => null() + ! + ! -- Initialize variables + inodata = 0 + nocdobj = 1 + allocate (this%ocdobj(nocdobj)) + do i = 1, nocdobj + call ocd_cr(ocdobjptr) + select case (i) + case (1) + call ocdobjptr%init_dbl('BUDGET', nullvec, dis, 'PRINT LAST ', & + 'COLUMNS 10 WIDTH 11 DIGITS 4 GENERAL ', & + this%iout, dnodata) + !case (2) + ! if stage is ever a dependent variable, we can activate following line + !call ocdobjptr%init_dbl('STAGE', head, dis, 'PRINT LAST ', & + ! 'COLUMNS 10 WIDTH 11 DIGITS 4 GENERAL ', & + ! this%iout, dnodata) + end select + this%ocdobj(i) = ocdobjptr + deallocate (ocdobjptr) + end do + ! + ! -- Read options or set defaults if this package not on + if (this%inunit > 0) then + call this%read_options() + end if + ! + ! -- Return + return + end subroutine oc_ar + +end module SnfOcModule diff --git a/src/SimulationCreate.f90 b/src/SimulationCreate.f90 index f5d70a1a057..8dc9e1a7136 100644 --- a/src/SimulationCreate.f90 +++ b/src/SimulationCreate.f90 @@ -250,6 +250,7 @@ subroutine models_create() ! -- modules use GwfModule, only: gwf_cr use GwtModule, only: gwt_cr + use SnfModule, only: snf_cr use ConstantsModule, only: LENMODELNAME ! -- dummy ! -- local @@ -282,6 +283,11 @@ subroutine models_create() call add_model(im, 'GWT6', mname) call gwt_cr(fname, im, modelname(im)) call add_dist_model(im) + case ('SNF6') + call parser%GetString(fname) + call add_model(im, 'SNF6', mname) + call snf_cr(fname, im, modelname(im)) + call add_dist_model(im) case default write (errmsg, '(4x,a,a)') & '****ERROR. UNKNOWN SIMULATION MODEL: ', & @@ -391,6 +397,7 @@ subroutine solution_groups_create() use BaseModelModule, only: BaseModelType use BaseExchangeModule, only: BaseExchangeType use NumericalSolutionModule, only: solution_create + use ExplicitSolutionModule, only: solution_explicit_create ! -- dummy ! -- local type(SolutionGroupType), pointer :: sgp @@ -498,6 +505,43 @@ subroutine solution_groups_create() ! end do ! + case ('EMS6') + ! + ! -- Initialize and increment counters + isoln = isoln + 1 + isgpsoln = isgpsoln + 1 + ! + ! -- Create the solution, retrieve from the list, and add to sgp + call parser%GetString(fname) + call solution_explicit_create(fname, isoln) + sp => GetBaseSolutionFromList(basesolutionlist, isoln) + call sgp%add_solution(isoln, sp) + ! + ! -- Add all of the models that are listed on this line to + ! the current solution (sp) + do + ! + ! -- Set istart and istop to encompass model name. Exit this + ! loop if there are no more models. + call parser%GetStringCaps(mname) + if (mname == '') exit + ! + ! -- Find the model id, and then get model + mid = ifind(modelname, mname) + if (mid <= 0) then + write (errmsg, '(a,a)') 'Error. Invalid modelname: ', & + trim(mname) + call store_error(errmsg) + call parser%StoreErrorUnit() + end if + mp => GetBaseModelFromList(basemodellist, mid) + ! + ! -- Add the model to the solution + call sp%add_model(mp) + mp%idsoln = isoln + ! + end do + ! case default write (errmsg, '(4x,a,a)') & '****ERROR. UNKNOWN SOLUTIONGROUP ENTRY: ', & diff --git a/src/Solution/ExplicitSolution.f90 b/src/Solution/ExplicitSolution.f90 new file mode 100644 index 00000000000..38746d1f6f0 --- /dev/null +++ b/src/Solution/ExplicitSolution.f90 @@ -0,0 +1,395 @@ +! This is the explicit solution module. + +module ExplicitSolutionModule + use KindModule, only: I4B, DP ! kluge note: DP needed for ttsoln + use TimerModule, only: code_timer + use ConstantsModule, only: LENMEMPATH, LENSOLUTIONNAME, MVALIDATE, & + MNORMAL, LINELENGTH, DZERO ! kluge note: DZERO needed for ttsoln + use MemoryHelperModule, only: create_mem_path + use BaseModelModule, only: BaseModelType, & + AddBaseModelToList, & + GetBaseModelFromList + use BaseExchangeModule, only: BaseExchangeType + use BaseSolutionModule, only: BaseSolutionType, AddBaseSolutionToList + use ListModule, only: ListType + use ListsModule, only: basesolutionlist + use SimVariablesModule, only: iout, isim_mode + use BlockParserModule, only: BlockParserType + + implicit none + private + + public :: solution_explicit_create + public :: ExplicitSolutionType + + type, extends(BaseSolutionType) :: ExplicitSolutionType + character(len=LENMEMPATH) :: memoryPath !< the path for storing solution variables in the memory manager + type(ListType), pointer :: modellist !< list of models in solution + integer(I4B), pointer :: id !< solution number + integer(I4B), pointer :: iu !< input file unit + real(DP), pointer :: ttsoln !< timer - total solution time + integer(I4B), pointer :: icnvg => null() !< convergence flag + type(BlockParserType) :: parser !< block parser object + contains + procedure :: sln_df + procedure :: sln_ar + procedure :: sln_calculate_delt + procedure :: sln_ad + procedure :: sln_ot + procedure :: sln_ca + procedure :: sln_fp + procedure :: sln_da + procedure :: add_model + procedure :: add_exchange + procedure :: get_models + procedure :: get_exchanges + procedure :: save + + procedure, private :: allocate_scalars + + ! Expose these for use through the BMI/XMI: + procedure, public :: prepareSolve + procedure, public :: solve + procedure, public :: finalizeSolve + + end type ExplicitSolutionType + +contains + +!> @ brief Create a new solution +!! +!! Create a new solution using the data in filename, assign this new +!! solution an id number and store the solution in the basesolutionlist. +!! Also open the filename for later reading. +!! +!< +subroutine solution_explicit_create(filename, id) + ! -- modules + use SimVariablesModule, only: iout + use InputOutputModule, only: getunit, openfile + ! -- dummy variables + character(len=*),intent(in) :: filename !< solution input file name + integer(I4B),intent(in) :: id !< solution id + ! -- local variables + integer(I4B) :: inunit + type(ExplicitSolutionType), pointer :: solution => null() + class(BaseSolutionType), pointer :: solbase => null() + character(len=LENSOLUTIONNAME) :: solutionname + ! + ! -- Create a new solution and add it to the basesolutionlist container + allocate(solution) + solbase => solution + write(solutionname,'(a, i0)') 'SLN_', id + ! + solution%name = solutionname + solution%memoryPath = create_mem_path(solutionname) + allocate(solution%modellist) + ! + call solution%allocate_scalars() + ! + call AddBaseSolutionToList(basesolutionlist, solbase) + ! + solution%id = id + ! + ! -- Open solution input file for reading later after problem size is known + ! Check to see if the file is already opened, which can happen when + ! running in single model mode + inquire(file=filename, number=inunit) + + if(inunit < 0) inunit = getunit() + solution%iu = inunit + write(iout,'(/a,a/)') ' Creating explicit solution (EMS): ', solution%name + call openfile(solution%iu, iout, filename, 'IMS') + ! + ! -- Initialize block parser + call solution%parser%Initialize(solution%iu, iout) + ! + ! -- return + return + end subroutine solution_explicit_create + + !> @ brief Allocate scalars + !! + !! Allocate scalars for a new solution. + !! + !< + subroutine allocate_scalars(this) + ! -- modules + use MemoryManagerModule, only: mem_allocate + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! + ! -- allocate scalars + call mem_allocate(this%id, 'ID', this%memoryPath) + call mem_allocate(this%iu, 'IU', this%memoryPath) + call mem_allocate(this%ttsoln, 'TTSOLN', this%memoryPath) + call mem_allocate(this%icnvg, 'ICNVG', this%memoryPath) + ! + ! -- initialize + this%id = 0 + this%iu = 0 + this%ttsoln = DZERO + this%icnvg = 0 + ! + ! -- return + return + end subroutine allocate_scalars + + subroutine sln_df(this) + class(ExplicitSolutionType) :: this + end subroutine + + subroutine sln_ar(this) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! + ! -- close ems input file + call this%parser%Clear() + ! + ! -- return + return + end subroutine sln_ar + + subroutine sln_calculate_delt(this) + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + end subroutine sln_calculate_delt + + !> @ brief Advance solution + !! + !! Advance solution. + !! + !< + subroutine sln_ad(this) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! + ! -- reset convergence flag + this%icnvg = 0 + + return + end subroutine sln_ad + + subroutine sln_ot(this) + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + end subroutine sln_ot + + subroutine sln_fp(this) + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + end subroutine sln_fp + + !> @ brief Deallocate solution + !! + !! Deallocate a solution. + !! + !< + subroutine sln_da(this) + ! -- modules + use MemoryManagerModule, only: mem_deallocate + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! + ! -- lists + call this%modellist%Clear() + deallocate(this%modellist) + ! + ! + ! -- Scalars + call mem_deallocate(this%id) + call mem_deallocate(this%iu) + call mem_deallocate(this%ttsoln) + call mem_deallocate(this%icnvg) + ! + ! -- return + return + end subroutine sln_da + + !> @ brief Solve solution + !! + !! Solve the models in this solution for kper and kstp. + !! + !< + subroutine sln_ca(this, isgcnvg, isuppress_output) + ! -- modules + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag + integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output + ! -- local variables + class(BaseModelType), pointer :: mp => null() + character(len=LINELENGTH) :: line + character(len=LINELENGTH) :: fmt + integer(I4B) :: im +! ------------------------------------------------------------------------------ + + ! advance the models and solution + call this%prepareSolve() + + select case (isim_mode) + case (MVALIDATE) + line = 'mode="validation" -- Skipping assembly and solution.' + fmt = "(/,1x,a,/)" + do im = 1, this%modellist%Count() + mp => GetBaseModelFromList(this%modellist, im) + call mp%model_message(line, fmt=fmt) + end do + case(MNORMAL) + + ! solve the models + call this%solve() + + ! finish up + call this%finalizeSolve(isgcnvg, isuppress_output) + end select + ! + ! -- return + return + + end subroutine sln_ca + + !> @ brief prepare to solve + !! + !! Prepare for the system solve by advancing the simulation. + !! + !< + subroutine prepareSolve(this) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! -- local variables + integer(I4B) :: im + class(BaseModelType), pointer :: mp => null() + + ! -- Model advance + do im = 1, this%modellist%Count() + mp => GetBaseModelFromList(this%modellist, im) + call mp%model_ad() + enddo + + ! advance solution + call this%sln_ad() + + end subroutine prepareSolve + + !> @ brief Build and solve the simulation + !! + !! Builds and solves the system for this explicit solution. + !! It roughly consists of the following steps + !! (1) ! kluge note: list steps + !! + !< + subroutine solve(this) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + ! -- local variables + class(BaseModelType), pointer :: mp => null() + integer(I4B) :: im + real(DP) :: ttsoln + ! + ! -- particle solve + call code_timer(0, ttsoln, this%ttsoln) ! kluge, needed? + do im = 1, this%modellist%Count() + mp => GetBaseModelFromList(this%modellist, im) + call mp%model_solve() + end do + call code_timer(1, ttsoln, this%ttsoln) ! kluge, needed? + ! + this%icnvg = 1 + + end subroutine solve + + !> @ brief finalize a solution + !! + !! Finalize the solution. Called after the outer iteration loop. + !! + !< + subroutine finalizeSolve(this, isgcnvg, isuppress_output) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag + integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output + ! -- local variables + integer(I4B) :: im + class(BaseModelType), pointer :: mp => null() + ! + ! -- Calculate flow for each model + do im=1,this%modellist%Count() + mp => GetBaseModelFromList(this%modellist, im) + call mp%model_cq(this%icnvg, isuppress_output) + enddo + ! + ! -- Budget terms for each model + do im=1,this%modellist%Count() + mp => GetBaseModelFromList(this%modellist, im) + call mp%model_bd(this%icnvg, isuppress_output) + enddo + ! + end subroutine finalizeSolve + + !> @ brief Save solution data to a file + !! + !! ! kluge note: description + !! + !< + subroutine save(this, filename) + ! -- modules + use InputOutputModule, only:getunit + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + character(len=*), intent(in) :: filename !< filename to save solution data + ! -- local variables + integer(I4B) :: inunit + ! + inunit = getunit() + open(unit=inunit,file=filename,status='unknown') + write(inunit,*) 'The save routine currently writes nothing' ! kluge + close(inunit) + ! + ! -- return + return + end subroutine save + + !> @ brief Add a model + !! + !! Add a model to this%modellist. + !! + !< + subroutine add_model(this, mp) + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + class(BaseModelType), pointer, intent(in) :: mp !< model instance + ! -- local variables + class(BaseModelType), pointer :: m => null() + ! + ! -- add a model + m => mp + call AddBaseModelToList(this%modellist, m) + ! + ! -- return + return + end subroutine add_model + + !> @brief Get a list of models + !! + !! Returns a pointer to the list of models in this solution. + !! + !< + function get_models(this) result(models) + ! -- return variable + type(ListType), pointer :: models !< pointer to the model list + ! -- dummy variables + class(ExplicitSolutionType) :: this !< ExplicitSolutionType instance + + models => this%modellist + + end function get_models + + subroutine add_exchange(this, exchange) + class(ExplicitSolutionType) :: this + class(BaseExchangeType), pointer, intent(in) :: exchange + end subroutine add_exchange + + function get_exchanges(this) result(exchanges) + class(ExplicitSolutionType) :: this + type(ListType), pointer :: exchanges + end function get_exchanges + +end module ExplicitSolutionModule \ No newline at end of file diff --git a/src/Utilities/Idm/InputDefinitionSelector.f90 b/src/Utilities/Idm/InputDefinitionSelector.f90 index c93d997db0d..3e0451e83dc 100644 --- a/src/Utilities/Idm/InputDefinitionSelector.f90 +++ b/src/Utilities/Idm/InputDefinitionSelector.f90 @@ -27,6 +27,12 @@ module InputDefinitionSelectorModule use GwtDspInputModule, only: gwt_dsp_param_definitions, & gwt_dsp_aggregate_definitions, & gwt_dsp_block_definitions + use SnfDislInputModule, only: snf_disl_param_definitions, & + snf_disl_aggregate_definitions, & + snf_disl_block_definitions + use SnfMmrInputModule, only: snf_mmr_param_definitions, & + snf_mmr_aggregate_definitions, & + snf_mmr_block_definitions implicit none private @@ -56,6 +62,10 @@ function param_definitions(component) result(input_definition) call set_pointer(input_definition, gwf_npf_param_definitions) case ('GWT/DSP') call set_pointer(input_definition, gwt_dsp_param_definitions) + case ('SNF/DISL') + call set_pointer(input_definition, snf_disl_param_definitions) + case ('SNF/MMR') + call set_pointer(input_definition, snf_mmr_param_definitions) case default write (warnmsg, '(a,a)') 'IDM Unsupported input type: ', trim(component) call store_warning(warnmsg) @@ -81,6 +91,10 @@ function aggregate_definitions(component) result(input_definition) call set_pointer(input_definition, gwf_npf_aggregate_definitions) case ('GWT/DSP') call set_pointer(input_definition, gwt_dsp_aggregate_definitions) + case ('SNF/DISL') + call set_pointer(input_definition, snf_disl_aggregate_definitions) + case ('SNF/MMR') + call set_pointer(input_definition, snf_mmr_aggregate_definitions) case default write (warnmsg, '(a,a)') 'IDM Unsupported input type: ', trim(component) call store_warning(warnmsg) @@ -106,6 +120,10 @@ function block_definitions(component) result(input_definition) call set_block_pointer(input_definition, gwf_npf_block_definitions) case ('GWT/DSP') call set_block_pointer(input_definition, gwt_dsp_block_definitions) + case ('SNF/DISL') + call set_block_pointer(input_definition, snf_disl_block_definitions) + case ('SNF/MMR') + call set_block_pointer(input_definition, snf_mmr_block_definitions) case default write (warnmsg, '(a,a)') 'IDM Unsupported input type: ', trim(component) call store_warning(warnmsg) diff --git a/src/Utilities/Idm/LoadMf6FileType.f90 b/src/Utilities/Idm/LoadMf6FileType.f90 index 24ff10d6e33..c5ea504ef27 100644 --- a/src/Utilities/Idm/LoadMf6FileType.f90 +++ b/src/Utilities/Idm/LoadMf6FileType.f90 @@ -116,9 +116,10 @@ subroutine parse_block(parser, mf6_input, iblock, mshape, iout) type(MemoryType), pointer :: mt ! ! -- disu vertices/cell2d blocks are contingent on NVERT dimension - if (mf6_input%file_type == 'DISU6') then + if (mf6_input%file_type == 'DISU6' .or. mf6_input%file_type == 'DISL6') then if (mf6_input%p_block_dfns(iblock)%blockname == 'VERTICES' .or. & - mf6_input%p_block_dfns(iblock)%blockname == 'CELL2D') then + mf6_input%p_block_dfns(iblock)%blockname == 'CELL2D' .or. & + mf6_input%p_block_dfns(iblock)%blockname == 'CELL1D') then call get_from_memorylist('NVERT', mf6_input%memoryPath, mt, found, & .false.) if (.not. found .or. mt%intsclr == 0) return @@ -741,7 +742,7 @@ subroutine set_model_shape(ftype, model_mempath, dis_mempath, model_shape) call mem_setptr(ndim1, 'NLAY', dis_mempath) call mem_setptr(ndim2, 'NCPL', dis_mempath) model_shape = [ndim1, ndim2] - case ('DISU6') + case ('DISU6', 'DISL6') call mem_allocate(model_shape, 1, 'MODEL_SHAPE', model_mempath) call mem_setptr(ndim1, 'NODES', dis_mempath) model_shape = [ndim1] diff --git a/src/meson.build b/src/meson.build index f50a97e7b7f..736c1d0a406 100644 --- a/src/meson.build +++ b/src/meson.build @@ -83,10 +83,18 @@ modflow_sources = files( 'Model' / 'GroundWaterTransport' / 'gwt1src1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1ssm1.f90', 'Model' / 'GroundWaterTransport' / 'gwt1uzt1.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1disl1.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1disl1idm.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1mmr1.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1mmr1idm.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1oc1.f90', + 'Model' / 'StreamNetworkFlow' / 'snf1flw1.f90', 'Model' / 'ModelUtilities' / 'BoundaryPackage.f90', 'Model' / 'ModelUtilities' / 'Connections.f90', 'Model' / 'ModelUtilities' / 'DiscretizationBase.f90', 'Model' / 'ModelUtilities' / 'DisvGeom.f90', + 'Model' / 'ModelUtilities' / 'DislGeom.f90', 'Model' / 'ModelUtilities' / 'GwfBuyInputData.f90', 'Model' / 'ModelUtilities' / 'GwfMvrPeriodData.f90', 'Model' / 'ModelUtilities' / 'GwfNpfOptions.f90', @@ -112,6 +120,7 @@ modflow_sources = files( 'Solution' / 'LinearMethods' / 'ims8misc.f90', 'Solution' / 'BaseSolution.f90', 'Solution' / 'NumericalSolution.f90', + 'Solution' / 'ExplicitSolution.f90', 'Solution' / 'SolutionGroup.f90', 'Timing' / 'ats.f90', 'Timing' / 'tdis.f90', diff --git a/utils/idmloader/scripts/dfn2f90.py b/utils/idmloader/scripts/dfn2f90.py index d0697f9eedd..10f00e688d0 100644 --- a/utils/idmloader/scripts/dfn2f90.py +++ b/utils/idmloader/scripts/dfn2f90.py @@ -448,6 +448,14 @@ def _source_file_footer(self, component, subcomponent): Path("../../../doc/mf6io/mf6ivar/dfn", "gwt-dsp.dfn"), Path("../../../src/Model/GroundWaterTransport", "gwt1dspidm.f90"), ], + [ + Path("../../../doc/mf6io/mf6ivar/dfn", "snf-disl.dfn"), + Path("../../../src/Model/StreamNetworkFlow", "snf1disl1idm.f90"), + ], + [ + Path("../../../doc/mf6io/mf6ivar/dfn", "snf-mmr.dfn"), + Path("../../../src/Model/StreamNetworkFlow", "snf1mmr1idm.f90"), + ], ] for dfn in dfns: