diff --git a/autotest/conftest.py b/autotest/conftest.py index f3bfadfa6a8..39633473390 100644 --- a/autotest/conftest.py +++ b/autotest/conftest.py @@ -128,18 +128,16 @@ def pytest_addoption(parser): def pytest_collection_modifyitems(config, items): - if config.getoption("--parallel"): - # --parallel given in cli: do not skip parallel tests - return - skip_parallel = pytest.mark.skip(reason="need --parallel option to run") - for item in items: - if "parallel" in item.keywords: - item.add_marker(skip_parallel) - - if config.getoption("--netcdf"): - # --netcdf given in cli: do not skip netcdf tests - return - skip_netcdf = pytest.mark.skip(reason="need --netcdf option to run") - for item in items: - if "netcdf" in item.keywords: - item.add_marker(skip_netcdf) + if not config.getoption("--parallel"): + skip_parallel = pytest.mark.skip( + reason="need --parallel option to run" + ) + for item in items: + if "parallel" in item.keywords: + item.add_marker(skip_parallel) + + if not config.getoption("--netcdf"): + skip_netcdf = pytest.mark.skip(reason="need --netcdf option to run") + for item in items: + if "netcdf" in item.keywords: + item.add_marker(skip_netcdf) diff --git a/autotest/test_netcdf_gwe_cnd.py b/autotest/test_netcdf_gwe_cnd.py index 4bb2f950396..598ff7cf515 100644 --- a/autotest/test_netcdf_gwe_cnd.py +++ b/autotest/test_netcdf_gwe_cnd.py @@ -20,15 +20,12 @@ msg += " pip install flopy" raise Exception(msg) -try: - import xarray as xa - import xugrid as xu -except ImportError: - pytest.skip("xarray and xugrid not found", allow_module_level=True) - from framework import TestFramework from test_gwe_cnd import cases +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") + def build_models(idx, test, export, gridded_input): from test_gwe_cnd import build_models as build @@ -162,7 +159,6 @@ def check_output(idx, test, export, gridded_input): nper = getattr(tdis, "nper").data nlay = getattr(dis, "nlay").data pd = getattr(tdis, "perioddata").array - print(pd) timestep = 0 for i in range(nper): for j in range(int(pd[i][1])): diff --git a/autotest/test_netcdf_gwf_disv.py b/autotest/test_netcdf_gwf_disv.py index bbbc5ee7989..ec716ddc73e 100644 --- a/autotest/test_netcdf_gwf_disv.py +++ b/autotest/test_netcdf_gwf_disv.py @@ -13,10 +13,8 @@ from framework import TestFramework from test_gwf_disv import cases -try: - import xugrid as xu -except ImportError: - pytest.skip("xugrid not found", allow_module_level=True) +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") wkt = ( 'PROJCS["NAD83 / UTM zone 18N", ' diff --git a/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py b/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py new file mode 100644 index 00000000000..d63d37fe3a9 --- /dev/null +++ b/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py @@ -0,0 +1,261 @@ +""" +NetCDF export test version of test_gwf_lak_wetlakbedarea02. This test compares +the temperature and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. +""" + +# Imports + +import os +import subprocess + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import TestFramework +from test_gwf_lak_wetlakbedarea02 import cases + +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") + + +def build_models(idx, test, export, gridded_input): + from test_gwf_lak_wetlakbedarea02 import build_models as build + + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = export + gwf.dis.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True + gwf.rch.export_array_netcdf = True + gwf.evt.export_array_netcdf = True + + name = cases[idx] + gwfname = "gwf-" + name + + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.dis, + filename=f"{gwfname}.dis.ncf", + ) + + return sim, dummy + + +def check_output(idx, test, export, gridded_input): + from test_gwf_lak_wetlakbedarea02 import check_output as check + + name = cases[idx] + gwfname = "gwf-" + name + + if gridded_input == "netcdf": + # re-run the simulation with model netcdf input + input_fname = f"{gwfname}.nc" + nc_fname = f"{gwfname}.{export}.nc" + subprocess.run( + ["mv", test.workspace / input_fname, test.workspace / nc_fname] + ) + + with open(test.workspace / f"{gwfname}.nam", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(" NEWTON\n") + f.write(f" EXPORT_NETCDF {export}\n") + f.write(f" NETCDF FILEIN {gwfname}.{export}.nc\n") + f.write("END options\n\n") + f.write("BEGIN packages\n") + f.write(f" DIS6 {gwfname}.dis dis\n") + f.write(f" NPF6 {gwfname}.npf npf\n") + f.write(f" STO6 {gwfname}.sto sto\n") + f.write(f" IC6 {gwfname}.ic ic\n") + f.write(f" CHD6 {gwfname}.chd chd_0\n") + f.write(f" RCH6 {gwfname}.rcha rcha_0\n") + f.write(f" EVT6 {gwfname}.evta evta_0\n") + f.write(f" LAK6 {gwfname}.lak lak-1\n") + f.write(f" OC6 {gwfname}.oc oc\n") + f.write("END packages\n") + + with open(test.workspace / f"{gwfname}.dis", "w") as f: + f.write("BEGIN options\n") + f.write(" LENGTH_UNITS feet\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write(f" NCF6 FILEIN {gwfname}.dis.ncf\n") + f.write("END options\n\n") + f.write("BEGIN dimensions\n") + f.write(" NLAY 6\n") + f.write(" NROW 17\n") + f.write(" NCOL 17\n") + f.write("END dimensions\n\n") + f.write("BEGIN griddata\n") + f.write(" delr NETCDF\n") + f.write(" delc NETCDF\n") + f.write(" top NETCDF\n") + f.write(" botm NETCDF\n") + f.write(" idomain NETCDF\n") + f.write("END griddata\n\n") + + with open(test.workspace / f"{gwfname}.ic", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" strt NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{gwfname}.npf", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_SPECIFIC_DISCHARGE\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" icelltype NETCDF\n") + f.write(" k NETCDF\n") + f.write(" k33 NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{gwfname}.rcha", "w") as f: + f.write("BEGIN options\n") + f.write(" READASARRAYS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" recharge NETCDF\n") + f.write("END period 1\n") + + with open(test.workspace / f"{gwfname}.evta", "w") as f: + f.write("BEGIN options\n") + f.write(" READASARRAYS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" surface NETCDF\n") + f.write(" rate NETCDF\n") + f.write(" depth NETCDF\n") + f.write("END period 1\n") + + success, buff = flopy.run_model( + test.targets["mf6"], + test.workspace / "mfsim.nam", + model_ws=test.workspace, + report=True, + ) + + assert success + test.success = success + + check(idx, test) + + # read transport results from GWF model + name = cases[idx] + gwfname = "gwf-" + name + + try: + # load heads + fname = gwfname + ".hds" + fpth = os.path.join(test.workspace, fname) + hobj = flopy.utils.HeadFile(fpth, precision="double") + heads = hobj.get_alldata() + except: + assert False, f'could not load headfile data from "{fpth}"' + + # Check NetCDF output + nc_fpth = os.path.join(test.workspace, f"{gwfname}.nc") + if export == "ugrid": + ds = xu.open_dataset(nc_fpth) + xds = ds.ugrid.to_dataset() + elif export == "structured": + xds = xa.open_dataset(nc_fpth) + + # Compare NetCDF head arrays with binary headfile temperatures + gwf = test.sims[0].gwf[0] + dis = getattr(gwf, "dis") + tdis = getattr(test.sims[0], "tdis") + nper = getattr(tdis, "nper").data + nlay = getattr(dis, "nlay").data + pd = getattr(tdis, "perioddata").array + timestep = 0 + for i in range(nper): + for j in range(int(pd[i][1])): + rec = hobj.get_data(kstpkper=(j, i)) + if export == "ugrid": + for l in range(nlay): + assert np.allclose( + np.array(rec[l]).flatten(), + # xds[f"head_l{l+1}"][timestep, :].data, + xds[f"head_l{l+1}"][timestep, :] + .fillna(1.00000000e30) + .data, + ), f"NetCDF-temperature comparison failure in timestep {timestep+1}" + timestep += 1 + elif export == "structured": + assert np.allclose( + # np.array(rec).flatten(), + np.array(rec), + xds["head"][timestep, :].fillna(1.00000000e30).data, + ), f"NetCDF-head comparison failure in timestep {timestep+1}" + timestep += 1 + + vlist = [ + "dis_delr", + "dis_delc", + "dis_top", + "dis_botm_l", + "dis_idomain_l", + "ic_strt_l", + "npf_icelltype_l", + "npf_k_l", + ] + + # Compare NetCDF package input arrays with FloPy arrays + gwf = test.sims[0].gwf[0] + for i, var in enumerate(vlist): + tokens = var.split("_", 1) + package_name = tokens[0] + array_name = tokens[1].split("_")[0] + package = getattr(gwf, package_name) + b = getattr(package, array_name).array + if export == "ugrid": + if var.endswith("_l"): + for l in range(nlay): + assert np.allclose( + np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data + ), f"NetCDF input array comparison failure, variable={var}{l+1}" + else: + assert np.allclose( + np.array(b).flatten(), xds[var].data + ), f"NetCDF input array comparison failure, variable={var}" + elif export == "structured": + var = var.replace("_l", "") + assert np.allclose( + # np.array(b).flatten(), xds[var].data + np.array(b), + xds[var].data, + ), f"NetCDF input array comparison failure, variable={var}" + + +@pytest.mark.netcdf +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +@pytest.mark.parametrize("export", ["ugrid", "structured"]) +@pytest.mark.parametrize("gridded_input", ["ascii", "netcdf"]) +def test_mf6model(idx, name, function_tmpdir, targets, export, gridded_input): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t, export, gridded_input), + check=lambda t: check_output(idx, t, export, gridded_input), + targets=targets, + ) + test.run() diff --git a/autotest/test_netcdf_gwf_rch01.py b/autotest/test_netcdf_gwf_rch01.py new file mode 100644 index 00000000000..78e1d5d60ea --- /dev/null +++ b/autotest/test_netcdf_gwf_rch01.py @@ -0,0 +1,258 @@ +""" +NetCDF export test version of test_gwf_rch01. This test compares +the temperature and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. +""" + +# Imports + +import os +import subprocess + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import TestFramework +from test_gwf_rch01 import cases + +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") + + +def build_models(idx, test, export, gridded_input): + from test_gwf_rch01 import build_models as build + + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = export + gwf.dis.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True + gwf.rch.export_array_netcdf = True + + name = "rch" + + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.dis, + filename=f"{name}.dis.ncf", + ) + + return sim, dummy + + +def check_output(idx, test, export, gridded_input): + from test_gwf_rch01 import check_output as check + + name = "rch" + + if gridded_input == "netcdf": + # re-run the simulation with model netcdf input + input_fname = f"{name}.nc" + nc_fname = f"{name}.{export}.nc" + subprocess.run( + ["mv", test.workspace / input_fname, test.workspace / nc_fname] + ) + + with open(test.workspace / f"{name}.nam", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(f" EXPORT_NETCDF {export}\n") + f.write(f" NETCDF FILEIN {name}.{export}.nc\n") + f.write("END options\n\n") + f.write("BEGIN packages\n") + f.write(f" DIS6 {name}.dis dis\n") + f.write(f" IC6 {name}.ic ic\n") + f.write(f" NPF6 {name}.npf npf\n") + f.write(f" STO6 {name}.sto sto\n") + f.write(f" CHD6 {name}.chd chd_0\n") + f.write(f" RCH6 {name}.rcha rcha_0\n") + f.write(f" OC6 {name}.oc oc\n") + f.write("END packages\n") + + with open(test.workspace / f"{name}.dis", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write(f" NCF6 FILEIN {name}.dis.ncf\n") + f.write("END options\n\n") + f.write("BEGIN dimensions\n") + f.write(" NLAY 2\n") + f.write(" NROW 1\n") + f.write(" NCOL 5\n") + f.write("END dimensions\n\n") + f.write("BEGIN griddata\n") + f.write(" delr NETCDF\n") + f.write(" delc NETCDF\n") + f.write(" top NETCDF\n") + f.write(" botm NETCDF\n") + f.write("END griddata\n\n") + + with open(test.workspace / f"{name}.ic", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" strt NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.npf", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" icelltype NETCDF\n") + f.write(" k NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.rcha", "w") as f: + f.write("BEGIN options\n") + f.write(" READASARRAYS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + if idx > 0: + f.write(" irch NETCDF\n") + f.write(" recharge NETCDF\n") + f.write("END period 1\n") + + success, buff = flopy.run_model( + test.targets["mf6"], + test.workspace / "mfsim.nam", + model_ws=test.workspace, + report=True, + ) + + assert success + test.success = success + + check(idx, test) + + # read transport results from GWF model + name = "rch" + + try: + # load heads + fpth = os.path.join(test.workspace, "rch.hds") + hobj = flopy.utils.HeadFile(fpth, precision="double") + heads = hobj.get_alldata() + except: + assert False, f'could not load headfile data from "{fpth}"' + + # Check NetCDF output + nc_fpth = os.path.join(test.workspace, f"{name}.nc") + if export == "ugrid": + ds = xu.open_dataset(nc_fpth) + xds = ds.ugrid.to_dataset() + elif export == "structured": + xds = xa.open_dataset(nc_fpth) + + # Compare NetCDF head arrays with binary headfile temperatures + gwf = test.sims[0].gwf[0] + dis = getattr(gwf, "dis") + tdis = getattr(test.sims[0], "tdis") + nper = getattr(tdis, "nper").data + nlay = getattr(dis, "nlay").data + pd = getattr(tdis, "perioddata").array + timestep = 0 + for i in range(nper): + for j in range(int(pd[i][1])): + rec = hobj.get_data(kstpkper=(j, i)) + if export == "ugrid": + for l in range(nlay): + assert np.allclose( + np.array(rec[l]).flatten(), + xds[f"head_l{l+1}"][timestep, :].data, + ), f"NetCDF-temperature comparison failure in timestep {timestep+1}" + timestep += 1 + elif export == "structured": + assert np.allclose( + # np.array(rec).flatten(), + np.array(rec), + xds["head"][timestep, :].data, + ), f"NetCDF-head comparison failure in timestep {timestep+1}" + timestep += 1 + + # compare recharge arrays + rch = getattr(gwf, "rcha_0") + irch = getattr(rch, "irch").array + recharge = getattr(rch, "recharge").array + if export == "ugrid": + rl1 = xds["rcha_0_recharge_l1_p1"].data.flatten() + rl2 = xds["rcha_0_recharge_l2_p1"].data.flatten() + elif export == "structured": + rl1 = xds["rcha_0_recharge_p1"].data[0].flatten() + rl2 = xds["rcha_0_recharge_p1"].data[1].flatten() + if idx == 1: + assert np.allclose( + np.array(irch).flatten() + 1, + xds["rcha_0_irch_p1"].data, + ), "NetCDF-irch comparison failure" + rarr = np.where(~np.isnan(rl1), rl1, rl2) + assert np.allclose( + np.array(recharge).flatten(), + rarr, + ), "NetCDF-recharge comparison failure" + + vlist = [ + "dis_delr", + "dis_delc", + "dis_top", + "dis_botm_l", + "ic_strt_l", + "npf_icelltype_l", + "npf_k_l", + ] + + # Compare NetCDF package input arrays with FloPy arrays + gwf = test.sims[0].gwf[0] + for i, var in enumerate(vlist): + tokens = var.split("_", 1) + package_name = tokens[0] + array_name = tokens[1].split("_")[0] + package = getattr(gwf, package_name) + b = getattr(package, array_name).array + if export == "ugrid": + if var.endswith("_l"): + for l in range(nlay): + assert np.allclose( + np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data + ), f"NetCDF input array comparison failure, variable={var}{l+1}" + else: + assert np.allclose( + np.array(b).flatten(), xds[var].data + ), f"NetCDF input array comparison failure, variable={var}" + elif export == "structured": + var = var.replace("_l", "") + assert np.allclose( + # np.array(b).flatten(), xds[var].data + np.array(b), + xds[var].data, + ), f"NetCDF input array comparison failure, variable={var}" + + +@pytest.mark.netcdf +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +@pytest.mark.parametrize("export", ["ugrid", "structured"]) +@pytest.mark.parametrize("gridded_input", ["ascii", "netcdf"]) +def test_mf6model(idx, name, function_tmpdir, targets, export, gridded_input): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t, export, gridded_input), + check=lambda t: check_output(idx, t, export, gridded_input), + targets=targets, + ) + test.run() diff --git a/autotest/test_netcdf_gwf_rch03.py b/autotest/test_netcdf_gwf_rch03.py new file mode 100644 index 00000000000..fc630d9582d --- /dev/null +++ b/autotest/test_netcdf_gwf_rch03.py @@ -0,0 +1,263 @@ +""" +NetCDF export test version of test_gwf_rch03. This test compares +the temperature and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. +""" + +# Imports + +import os +import subprocess + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import TestFramework +from test_gwf_rch03 import cases + +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") + + +def build_models(idx, test, export, gridded_input): + from test_gwf_rch03 import build_models as build + + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = export + gwf.dis.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True + gwf.rch.export_array_netcdf = True + print(gwf.rch) + + # name = "gwf-" + cases[idx] + name = "rch" + + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.dis, + filename=f"{name}.dis.ncf", + ) + + return sim, dummy + + +def check_output(idx, test, export, gridded_input): + from test_gwf_rch03 import check_output as check + + name = "rch" + + if gridded_input == "netcdf": + # re-run the simulation with model netcdf input + input_fname = f"{name}.nc" + nc_fname = f"{name}.{export}.nc" + subprocess.run( + ["mv", test.workspace / input_fname, test.workspace / nc_fname] + ) + + with open(test.workspace / f"{name}.nam", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(f" EXPORT_NETCDF {export}\n") + f.write(f" NETCDF FILEIN {name}.{export}.nc\n") + f.write("END options\n\n") + f.write("BEGIN packages\n") + f.write(f" DIS6 {name}.dis dis\n") + f.write(f" IC6 {name}.ic ic\n") + f.write(f" NPF6 {name}.npf npf\n") + f.write(f" CHD6 {name}.chd chd_0\n") + f.write(f" RCH6 {name}.rcha rcha_0\n") + f.write(f" OC6 {name}.oc oc\n") + f.write("END packages\n") + + with open(test.workspace / f"{name}.dis", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write(f" NCF6 FILEIN {name}.dis.ncf\n") + f.write("END options\n\n") + f.write("BEGIN dimensions\n") + f.write(" NLAY 2\n") + f.write(" NROW 4\n") + f.write(" NCOL 5\n") + f.write("END dimensions\n\n") + f.write("BEGIN griddata\n") + f.write(" delr NETCDF\n") + f.write(" delc NETCDF\n") + f.write(" top NETCDF\n") + f.write(" botm NETCDF\n") + f.write(" idomain NETCDF\n") + f.write("END griddata\n\n") + + with open(test.workspace / f"{name}.ic", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" strt NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.npf", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" icelltype NETCDF\n") + f.write(" k NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.rcha", "w") as f: + f.write("BEGIN options\n") + f.write(" READASARRAYS\n") + f.write(" PRINT_FLOWS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" irch NETCDF\n") + f.write(" recharge NETCDF\n") + f.write("END period 1\n") + + success, buff = flopy.run_model( + test.targets["mf6"], + test.workspace / "mfsim.nam", + model_ws=test.workspace, + report=True, + ) + + assert success + test.success = success + + check(idx, test) + + # read transport results from GWF model + name = "rch" + + try: + # load heads + fpth = os.path.join(test.workspace, "rch.hds") + hobj = flopy.utils.HeadFile(fpth, precision="double") + heads = hobj.get_alldata() + except: + assert False, f'could not load headfile data from "{fpth}"' + + # Check NetCDF output + nc_fpth = os.path.join(test.workspace, f"{name}.nc") + if export == "ugrid": + ds = xu.open_dataset(nc_fpth) + xds = ds.ugrid.to_dataset() + elif export == "structured": + xds = xa.open_dataset(nc_fpth) + + # Compare NetCDF head arrays with binary headfile temperatures + gwf = test.sims[0].gwf[0] + dis = getattr(gwf, "dis") + tdis = getattr(test.sims[0], "tdis") + nper = getattr(tdis, "nper").data + nlay = getattr(dis, "nlay").data + pd = getattr(tdis, "perioddata").array + print(pd) + timestep = 0 + for i in range(nper): + for j in range(int(pd[i][1])): + rec = hobj.get_data(kstpkper=(j, i)) + if export == "ugrid": + for l in range(nlay): + assert np.allclose( + np.array(rec[l]).flatten(), + xds[f"head_l{l+1}"][timestep, :] + .fillna(1.00000000e30) + .data, + ), f"NetCDF-temperature comparison failure in timestep {timestep+1}" + timestep += 1 + elif export == "structured": + assert np.allclose( + # np.array(rec).flatten(), + np.array(rec), + xds["head"][timestep, :].fillna(1.00000000e30).data, + ), f"NetCDF-head comparison failure in timestep {timestep+1}" + timestep += 1 + + # compare recharge arrays + rch = getattr(gwf, "rcha_0") + irch = getattr(rch, "irch").array + recharge = getattr(rch, "recharge").array + if export == "ugrid": + rl1 = xds["rcha_0_recharge_l1_p1"].data.flatten() + rl2 = xds["rcha_0_recharge_l2_p1"].data.flatten() + elif export == "structured": + rl1 = xds["rcha_0_recharge_p1"].data[0].flatten() + rl2 = xds["rcha_0_recharge_p1"].data[1].flatten() + assert np.allclose( + np.array(irch).flatten() + 1, + xds["rcha_0_irch_p1"].data, + ), "NetCDF-irch comparison failure" + rarr = np.where(~np.isnan(rl1), rl1, rl2) + assert np.allclose( + np.array(recharge).flatten(), + rarr, + ), "NetCDF-recharge comparison failure" + + vlist = [ + "dis_delr", + "dis_delc", + "dis_top", + "dis_botm_l", + "dis_idomain_l", + "ic_strt_l", + "npf_icelltype_l", + "npf_k_l", + ] + + # Compare NetCDF package input arrays with FloPy arrays + gwf = test.sims[0].gwf[0] + for i, var in enumerate(vlist): + tokens = var.split("_", 1) + package_name = tokens[0] + array_name = tokens[1].split("_")[0] + package = getattr(gwf, package_name) + b = getattr(package, array_name).array + if export == "ugrid": + if var.endswith("_l"): + for l in range(nlay): + assert np.allclose( + np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data + ), f"NetCDF input array comparison failure, variable={var}{l+1}" + else: + assert np.allclose( + np.array(b).flatten(), xds[var].data + ), f"NetCDF input array comparison failure, variable={var}" + elif export == "structured": + var = var.replace("_l", "") + assert np.allclose( + # np.array(b).flatten(), xds[var].data + np.array(b), + xds[var].data, + ), f"NetCDF input array comparison failure, variable={var}" + + +@pytest.mark.netcdf +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +@pytest.mark.parametrize("export", ["ugrid", "structured"]) +@pytest.mark.parametrize("gridded_input", ["ascii", "netcdf"]) +def test_mf6model(idx, name, function_tmpdir, targets, export, gridded_input): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t, export, gridded_input), + check=lambda t: check_output(idx, t, export, gridded_input), + targets=targets, + ) + test.run() diff --git a/autotest/test_netcdf_gwf_sto01.py b/autotest/test_netcdf_gwf_sto01.py index 8de196a3eef..da9744b3ec4 100644 --- a/autotest/test_netcdf_gwf_sto01.py +++ b/autotest/test_netcdf_gwf_sto01.py @@ -11,11 +11,8 @@ from framework import TestFramework from test_gwf_sto01 import cases -try: - import xarray as xa - import xugrid as xu -except ImportError: - pytest.skip("xarray and xugrid not found", allow_module_level=True) +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") htol = [None for _ in range(len(cases))] diff --git a/autotest/test_netcdf_gwf_vsc03_sfr.py b/autotest/test_netcdf_gwf_vsc03_sfr.py new file mode 100644 index 00000000000..2e6602bb4cf --- /dev/null +++ b/autotest/test_netcdf_gwf_vsc03_sfr.py @@ -0,0 +1,293 @@ +""" +NetCDF export test version of test_gwf_vsc03_sfr. This test compares +the temperature and input arrays in the the NetCDF file to those +in the FloPy binary output head file and package data objects. +""" + +# Imports + +import os +import subprocess + +import numpy as np +import pytest + +try: + import flopy +except: + msg = "Error. FloPy package is not available.\n" + msg += "Try installing using the following command:\n" + msg += " pip install flopy" + raise Exception(msg) + +from framework import TestFramework +from test_gwf_vsc03_sfr import cases + +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") + + +def build_models(idx, test, export, gridded_input): + from test_gwf_vsc03_sfr import build_models as build + + sim, dummy = build(idx, test) + sim.tdis.start_date_time = "2041-01-01T00:00:00-05:00" + gwf = sim.gwf[0] + gwf.name_file.export_netcdf = export + gwf.dis.export_array_netcdf = True + gwf.ic.export_array_netcdf = True + gwf.npf.export_array_netcdf = True + gwf.sto.export_array_netcdf = True + gwf.rcha.export_array_netcdf = True + assert gwf.rcha + assert gwf.sto + + name = "gwf-" + cases[idx] + + # netcdf config + ncf = flopy.mf6.ModflowUtlncf( + gwf.dis, + filename=f"{name}.dis.ncf", + ) + + return sim, dummy + + +def check_output(idx, test, export, gridded_input): + from test_gwf_vsc03_sfr import check_output as check + + name = "gwf-" + test.name + + if gridded_input == "netcdf": + # re-run the simulation with model netcdf input + input_fname = f"{name}.nc" + nc_fname = f"{name}.{export}.nc" + subprocess.run( + ["mv", test.workspace / input_fname, test.workspace / nc_fname] + ) + + with open(test.workspace / f"{name}.nam", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_FLOWS\n") + f.write(" NEWTON\n") + f.write(f" EXPORT_NETCDF {export}\n") + f.write(f" NETCDF FILEIN {name}.{export}.nc\n") + f.write("END options\n\n") + f.write("BEGIN packages\n") + f.write(f" DIS6 {name}.dis dis\n") + f.write(f" NPF6 {name}.npf npf\n") + f.write(f" STO6 {name}.sto sto\n") + f.write(f" IC6 {name}.ic ic\n") + if idx == 1: + f.write(" VSC6 gwf-vsc-sfr01.vsc vsc\n") + f.write(f" OC6 {name}.oc oc\n") + f.write(f" RCH6 {name}.rcha rcha-1\n") + f.write(f" EVT6 {name}.evt evt-1\n") + f.write(f" SFR6 {name}.sfr sfr-1\n") + f.write("END packages\n") + + with open(test.workspace / f"{name}.dis", "w") as f: + f.write("BEGIN options\n") + f.write(" LENGTH_UNITS m\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write(f" NCF6 FILEIN {name}.dis.ncf\n") + f.write("END options\n\n") + f.write("BEGIN dimensions\n") + f.write(" NLAY 1\n") + f.write(" NROW 60\n") + f.write(" NCOL 200\n") + f.write("END dimensions\n\n") + f.write("BEGIN griddata\n") + f.write(" delr NETCDF\n") + f.write(" delc NETCDF\n") + f.write(" top NETCDF\n") + f.write(" botm NETCDF\n") + f.write("END griddata\n\n") + + with open(test.workspace / f"{name}.ic", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" strt NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.npf", "w") as f: + f.write("BEGIN options\n") + f.write(" SAVE_SPECIFIC_DISCHARGE\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" icelltype NETCDF\n") + f.write(" k NETCDF\n") + f.write("END griddata\n") + + with open(test.workspace / f"{name}.sto", "w") as f: + f.write("BEGIN options\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN griddata\n") + f.write(" iconvert NETCDF\n") + f.write(" ss NETCDF\n") + f.write(" sy NETCDF\n") + f.write("END griddata\n") + f.write("BEGIN period 1\n") + f.write(" STEADY-STATE\n") + f.write("END period 1\n") + f.write("BEGIN period 2\n") + f.write(" TRANSIENT\n") + f.write("END period 2\n") + + with open(test.workspace / f"{name}.rcha", "w") as f: + f.write("BEGIN options\n") + f.write(" READASARRAYS\n") + f.write(" auxiliary TEMPERATURE\n") + f.write(" PRINT_FLOWS\n") + f.write(" EXPORT_ARRAY_NETCDF\n") + f.write("END options\n\n") + f.write("BEGIN period 1\n") + f.write(" irch NETCDF\n") + f.write(" recharge NETCDF\n") + f.write(" TEMPERATURE NETCDF\n") + f.write("END period 1\n") + + success, buff = flopy.run_model( + test.targets["mf6"], + test.workspace / "mfsim.nam", + model_ws=test.workspace, + report=True, + ) + + assert success + test.success = success + + check(idx, test) + + # read transport results from GWF model + name = "gwf-" + cases[idx] + + try: + # load heads + fpth = os.path.join(test.workspace, f"{name}.hds") + hobj = flopy.utils.HeadFile(fpth, precision="double") + heads = hobj.get_alldata() + except: + assert False, f'could not load headfile data from "{fpth}"' + + # Check NetCDF output + nc_fpth = os.path.join(test.workspace, f"{name}.nc") + if export == "ugrid": + ds = xu.open_dataset(nc_fpth) + xds = ds.ugrid.to_dataset() + elif export == "structured": + xds = xa.open_dataset(nc_fpth) + + # Compare NetCDF head arrays with binary headfile temperatures + gwf = test.sims[0].gwf[0] + dis = getattr(gwf, "dis") + tdis = getattr(test.sims[0], "tdis") + nper = getattr(tdis, "nper").data + nlay = getattr(dis, "nlay").data + pd = getattr(tdis, "perioddata").array + timestep = 0 + for i in range(nper): + for j in range(int(pd[i][1])): + rec = hobj.get_data(kstpkper=(j, i)) + if export == "ugrid": + for l in range(nlay): + assert np.allclose( + np.array(rec[l]).flatten(), + # xds[f"head_l{l+1}"][timestep, :].data, + xds[f"head_l{l+1}"][timestep, :] + .fillna(1.00000000e30) + .data, + ), f"NetCDF-temperature comparison failure in timestep {timestep+1}" + timestep += 1 + elif export == "structured": + assert np.allclose( + # np.array(rec).flatten(), + np.array(rec), + xds["head"][timestep, :].fillna(1.00000000e30).data, + ), f"NetCDF-head comparison failure in timestep {timestep+1}" + timestep += 1 + + # compare recharge arrays + rch = getattr(gwf, "rcha-1") + irch = getattr(rch, "irch").array + recharge = getattr(rch, "recharge").array + aux = getattr(rch, "aux").array + if export == "ugrid": + rarr = xds["rcha-1_recharge_l1_p1"].data.flatten() + auxarr = xds["rcha-1_auxvar_l1_p1a1"].data.flatten() + elif export == "structured": + rarr = xds["rcha-1_recharge_p1"].data[0].flatten() + auxarr = xds["rcha-1_auxvar_p1a1"].data[0].flatten() + assert np.allclose( + np.array(irch[0]).flatten() + 1, + xds["rcha-1_irch_p1"].data, + ), "NetCDF-irch comparison failure" + assert np.allclose( + np.array(recharge[0]).flatten(), + rarr, + ), "NetCDF-recharge comparison failure" + assert np.allclose( + np.array(aux[0]).flatten(), + auxarr, + ), "NetCDF-auxvar comparison failure" + + vlist = [ + "dis_delr", + "dis_delc", + "dis_top", + "dis_botm_l", + "ic_strt_l", + "npf_icelltype_l", + "npf_k_l", + "sto_iconvert_l", + "sto_ss_l", + "sto_sy_l", + ] + + # Compare NetCDF package input arrays with FloPy arrays + gwf = test.sims[0].gwf[0] + for i, var in enumerate(vlist): + tokens = var.split("_", 1) + package_name = tokens[0] + array_name = tokens[1].split("_")[0] + package = getattr(gwf, package_name) + b = getattr(package, array_name).array + if export == "ugrid": + if var.endswith("_l"): + for l in range(nlay): + assert np.allclose( + np.array(b[l]).flatten(), xds[f"{var}{l+1}"].data + ), f"NetCDF input array comparison failure, variable={var}{l+1}" + else: + assert np.allclose( + np.array(b).flatten(), xds[var].data + ), f"NetCDF input array comparison failure, variable={var}" + elif export == "structured": + var = var.replace("_l", "") + assert np.allclose( + # np.array(b).flatten(), xds[var].data + np.array(b), + xds[var].data, + ), f"NetCDF input array comparison failure, variable={var}" + + +@pytest.mark.netcdf +@pytest.mark.parametrize( + "idx, name", + list(enumerate(cases)), +) +@pytest.mark.parametrize("export", ["ugrid", "structured"]) +@pytest.mark.parametrize("gridded_input", ["ascii", "netcdf"]) +def test_mf6model(idx, name, function_tmpdir, targets, export, gridded_input): + test = TestFramework( + name=name, + workspace=function_tmpdir, + build=lambda t: build_models(idx, t, export, gridded_input), + check=lambda t: check_output(idx, t, export, gridded_input), + targets=targets, + ) + test.run() diff --git a/autotest/test_netcdf_gwt_dsp01.py b/autotest/test_netcdf_gwt_dsp01.py index c1ef491d637..7360a7956bc 100644 --- a/autotest/test_netcdf_gwt_dsp01.py +++ b/autotest/test_netcdf_gwt_dsp01.py @@ -11,11 +11,8 @@ from framework import TestFramework from test_gwt_dsp01 import cases, xt3d -try: - import xarray as xa - import xugrid as xu -except ImportError: - pytest.skip("xarray and xugrid not found", allow_module_level=True) +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") def build_models(idx, test, export, gridded_input): diff --git a/autotest/test_netcdf_gwt_prudic2004t2.py b/autotest/test_netcdf_gwt_prudic2004t2.py index fc9cfee8450..24efe36ac5a 100644 --- a/autotest/test_netcdf_gwt_prudic2004t2.py +++ b/autotest/test_netcdf_gwt_prudic2004t2.py @@ -11,11 +11,8 @@ from framework import TestFramework from test_gwt_prudic2004t2 import cases -try: - import xarray as xa - import xugrid as xu -except ImportError: - pytest.skip("xarray and xugrid not found", allow_module_level=True) +xa = pytest.importorskip("xarray") +xu = pytest.importorskip("xugrid") def build_models(idx, test, export, gridded_input): diff --git a/doc/mf6io/mf6ivar/dfn/gwf-evta.dfn b/doc/mf6io/mf6ivar/dfn/gwf-evta.dfn index 2ee518627ee..ab4a44c276f 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-evta.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-evta.dfn @@ -141,6 +141,14 @@ optional false longname obs6 input filename description REPLACE obs6_filename {'{#1}': 'Evapotranspiration'} +block options +name export_array_netcdf +type keyword +reader urword +optional true +mf6internal export_nc +longname export array variables to netcdf output files. +description keyword that specifies input griddata arrays should be written to the model output netcdf file. # --------------------- gwf evta period --------------------- @@ -201,6 +209,6 @@ type double precision shape (ncol*nrow; ncpl) reader readarray time_series true -longname auxiliary variable iaux +longname evapotranspiration auxiliary variable iaux description is an array of values for auxiliary variable AUX(IAUX), where iaux is a value from 1 to NAUX, and AUX(IAUX) must be listed as part of the auxiliary variables. A separate array can be specified for each auxiliary variable. If an array is not specified for an auxiliary variable, then a value of zero is assigned. If the value specified here for the auxiliary variable is the same as auxmultname, then the evapotranspiration rate will be multiplied by this array. mf6internal auxvar diff --git a/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn b/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn index f25fbfbd320..1441a0988b9 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-rcha.dfn @@ -141,6 +141,14 @@ optional false longname obs6 input filename description REPLACE obs6_filename {'{#1}': 'Recharge'} +block options +name export_array_netcdf +type keyword +reader urword +optional true +mf6internal export_nc +longname export array variables to netcdf output files. +description keyword that specifies input griddata arrays should be written to the model output netcdf file. # --------------------- gwf rcha period --------------------- @@ -184,6 +192,6 @@ shape (ncol*nrow; ncpl) reader readarray time_series true optional true -longname auxiliary variable iaux +longname recharge auxiliary variable iaux description is an array of values for auxiliary variable aux(iaux), where iaux is a value from 1 to naux, and aux(iaux) must be listed as part of the auxiliary variables. A separate array can be specified for each auxiliary variable. If an array is not specified for an auxiliary variable, then a value of zero is assigned. If the value specified here for the auxiliary variable is the same as auxmultname, then the recharge array will be multiplied by this array. mf6internal auxvar diff --git a/src/Idm/gwf-evtaidm.f90 b/src/Idm/gwf-evtaidm.f90 index 4820a429347..d3e74339ca8 100644 --- a/src/Idm/gwf-evtaidm.f90 +++ b/src/Idm/gwf-evtaidm.f90 @@ -26,6 +26,7 @@ module GwfEvtaInputModule logical :: obs_filerecord = .false. logical :: obs6 = .false. logical :: obs6_filename = .false. + logical :: export_nc = .false. logical :: ievt = .false. logical :: surface = .false. logical :: rate = .false. @@ -293,6 +294,24 @@ module GwfEvtaInputModule .false. & ! timeseries ) + type(InputParamDefinitionType), parameter :: & + gwfevta_export_nc = InputParamDefinitionType & + ( & + 'GWF', & ! component + 'EVTA', & ! subcomponent + 'OPTIONS', & ! block + 'EXPORT_ARRAY_NETCDF', & ! tag name + 'EXPORT_NC', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'export array variables to netcdf output files.', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + type(InputParamDefinitionType), parameter :: & gwfevta_ievt = InputParamDefinitionType & ( & @@ -375,7 +394,7 @@ module GwfEvtaInputModule 'AUXVAR', & ! fortran variable 'DOUBLE2D', & ! type 'NAUX NCPL', & ! shape - 'auxiliary variable iaux', & ! longname + 'evapotranspiration auxiliary variable iaux', & ! longname .true., & ! required .false., & ! multi-record .false., & ! preserve case @@ -400,6 +419,7 @@ module GwfEvtaInputModule gwfevta_obs_filerecord, & gwfevta_obs6, & gwfevta_obs6_filename, & + gwfevta_export_nc, & gwfevta_ievt, & gwfevta_surface, & gwfevta_rate, & diff --git a/src/Idm/gwf-rchaidm.f90 b/src/Idm/gwf-rchaidm.f90 index a80daa67182..009149a6cf3 100644 --- a/src/Idm/gwf-rchaidm.f90 +++ b/src/Idm/gwf-rchaidm.f90 @@ -26,6 +26,7 @@ module GwfRchaInputModule logical :: obs_filerecord = .false. logical :: obs6 = .false. logical :: obs6_filename = .false. + logical :: export_nc = .false. logical :: irch = .false. logical :: recharge = .false. logical :: auxvar = .false. @@ -291,6 +292,24 @@ module GwfRchaInputModule .false. & ! timeseries ) + type(InputParamDefinitionType), parameter :: & + gwfrcha_export_nc = InputParamDefinitionType & + ( & + 'GWF', & ! component + 'RCHA', & ! subcomponent + 'OPTIONS', & ! block + 'EXPORT_ARRAY_NETCDF', & ! tag name + 'EXPORT_NC', & ! fortran variable + 'KEYWORD', & ! type + '', & ! shape + 'export array variables to netcdf output files.', & ! longname + .false., & ! required + .false., & ! multi-record + .false., & ! preserve case + .false., & ! layered + .false. & ! timeseries + ) + type(InputParamDefinitionType), parameter :: & gwfrcha_irch = InputParamDefinitionType & ( & @@ -337,7 +356,7 @@ module GwfRchaInputModule 'AUXVAR', & ! fortran variable 'DOUBLE2D', & ! type 'NAUX NCPL', & ! shape - 'auxiliary variable iaux', & ! longname + 'recharge auxiliary variable iaux', & ! longname .false., & ! required .false., & ! multi-record .false., & ! preserve case @@ -362,6 +381,7 @@ module GwfRchaInputModule gwfrcha_obs_filerecord, & gwfrcha_obs6, & gwfrcha_obs6_filename, & + gwfrcha_export_nc, & gwfrcha_irch, & gwfrcha_recharge, & gwfrcha_auxvar & diff --git a/src/Utilities/Export/DisNCMesh.f90 b/src/Utilities/Export/DisNCMesh.f90 index 05ebbe18d80..8dfd56f155d 100644 --- a/src/Utilities/Export/DisNCMesh.f90 +++ b/src/Utilities/Export/DisNCMesh.f90 @@ -10,10 +10,15 @@ module MeshDisModelModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENBIGLINE, LENCOMPONENTNAME, & LENMEMPATH + use SimVariablesModule, only: errmsg + use SimModule, only: store_error, store_error_filename use MemoryManagerModule, only: mem_setptr use InputDefinitionModule, only: InputParamDefinitionType use CharacterStringModule, only: CharacterStringType - use MeshModelModule, only: Mesh2dModelType, MeshNCDimIdType, MeshNCVarIdType + use MeshModelModule, only: Mesh2dModelType, MeshNCDimIdType, MeshNCVarIdType, & + ncvar_chunk, ncvar_deflate, ncvar_gridmap, & + ncvar_mf6attr, export_varname + use NCModelExportModule, only: export_longname use DisModule, only: DisType use NetCDFCommonModule, only: nf_verify use netcdf @@ -33,6 +38,9 @@ module MeshDisModelModule procedure :: df procedure :: step procedure :: export_input_array + procedure :: package_step_ilayer + procedure :: package_step + procedure :: export_layer_3d procedure :: define_dim procedure :: add_mesh_data end type Mesh2dDisExportType @@ -172,11 +180,177 @@ subroutine step(this) nullify (dbl2d) end subroutine step + !> @brief netcdf export package dynamic input with ilayer index variable + !< + subroutine package_step_ilayer(this, export_pkg, ilayer_varname, ilayer) + use ConstantsModule, only: DNODATA, DZERO + use TdisModule, only: kper + use DefinitionSelectModule, only: get_param_definition_type + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + character(len=*), intent(in) :: ilayer_varname + integer(I4B), intent(in) :: ilayer + ! -- local + type(InputParamDefinitionType), pointer :: idt + integer(I4B), dimension(:), pointer, contiguous :: int1d + real(DP), dimension(:), pointer, contiguous :: dbl1d + real(DP), dimension(:, :), pointer, contiguous :: dbl2d + integer(I4B), dimension(:), pointer, contiguous :: ialayer + real(DP), dimension(:), contiguous, pointer :: dbl1d_ptr + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: n, iparam, nvals + logical(LGP) :: ilayer_read + ! + ! -- initialize + nullify (ialayer) + ilayer_read = .false. + ! + ! -- set pointer to ilayer variable + call mem_setptr(ialayer, export_pkg%param_names(ilayer), & + export_pkg%mf6_input%mempath) + ! + ! -- check if layer index variable was read + if (export_pkg%param_reads(ilayer)%invar == 1) then + ilayer_read = .true. + end if + ! + ! -- export defined period input + do iparam = 1, export_pkg%nparam + ! + ! -- check if variable was read this period + if (export_pkg%param_reads(iparam)%invar < 1) cycle + ! + ! -- set input definition + idt => & + get_param_definition_type(export_pkg%mf6_input%param_dfns, & + export_pkg%mf6_input%component_type, & + export_pkg%mf6_input%subcomponent_type, & + 'PERIOD', export_pkg%param_names(iparam), '') + ! + ! -- set variable name and input string + nc_varname = trim(export_pkg%mf6_input%subcomponent_name)//'_'// & + trim(idt%mf6varname) + input_attr = this%input_attribute(export_pkg%mf6_input%subcomponent_name, & + idt) + ! + ! -- export arrays + select case (idt%datatype) + case ('INTEGER1D') + call mem_setptr(int1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call nc_export_int1d(this%ncid, this%dim_ids, this%x_dim, this%y_dim, & + this%var_ids, this%dis, int1d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, & + idt%tagname, this%gridmap_name, idt%shape, & + idt%longname, input_attr, this%deflate, & + this%shuffle, this%chunk_face, kper, this%nc_fname) + case ('DOUBLE1D') + call mem_setptr(dbl1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call this%export_layer_3d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr) + case ('DOUBLE2D') + call mem_setptr(dbl2d, idt%mf6varname, export_pkg%mf6_input%mempath) + nvals = this%dis%ncol * this%dis%nrow + ! + do n = 1, size(dbl2d, dim=1) !naux + dbl1d_ptr(1:nvals) => dbl2d(n, :) + if (all(dbl1d_ptr == DZERO)) then + else + call this%export_layer_3d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d_ptr, nc_varname, input_attr, n) + end if + end do + case default + ! + errmsg = 'EXPORT ilayaer unsupported datatype='//trim(idt%datatype) + call store_error(errmsg, .true.) + end select + end do + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + ! + ! -- return + return + end subroutine package_step_ilayer + + !> @brief netcdf export package dynamic input + !< + subroutine package_step(this, export_pkg) + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + errmsg = 'NetCDF period export not supported for model='// & + trim(this%modelname)//', package='// & + trim(export_pkg%mf6_input%subcomponent_name) + call store_error(errmsg, .true.) + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + end subroutine package_step + + !> @brief export layer variable as full grid + !< + subroutine export_layer_3d(this, export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr, iaux) + use ConstantsModule, only: DNODATA, DZERO + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + type(InputParamDefinitionType), pointer, intent(in) :: idt + logical(LGP), intent(in) :: ilayer_read + integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ialayer + real(DP), dimension(:), pointer, contiguous, intent(in) :: dbl1d + character(len=*), intent(in) :: nc_varname + character(len=*), intent(in) :: input_attr + integer(I4B), optional, intent(in) :: iaux + ! -- local + real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d + integer(I4B) :: n, i, j, k, nvals, idxaux + real(DP), dimension(:, :), contiguous, pointer :: dbl2d_ptr + ! + ! -- initialize + idxaux = 0 + if (present(iaux)) then + idxaux = iaux + end if + + allocate (dbl3d(export_pkg%mshape(3), export_pkg%mshape(2), & + export_pkg%mshape(1))) + ! + if (ilayer_read) then + do k = 1, size(dbl3d, dim=3) + n = 0 + do i = 1, size(dbl3d, dim=2) + do j = 1, size(dbl3d, dim=1) + n = n + 1 + if (ialayer(n) == k) then + dbl3d(j, i, k) = dbl1d(n) + else + dbl3d(j, i, k) = DNODATA + end if + end do + end do + end do + else + dbl3d = DNODATA + nvals = export_pkg%mshape(3) * export_pkg%mshape(2) + dbl2d_ptr(1:export_pkg%mshape(3), 1:export_pkg%mshape(2)) => dbl1d(1:nvals) + dbl3d(:, :, 1) = dbl2d_ptr(:, :) + end if + ! + call nc_export_dbl3d(this%ncid, this%dim_ids, this%var_ids, this%dis, dbl3d, & + nc_varname, export_pkg%mf6_input%subcomponent_name, & + idt%tagname, this%gridmap_name, idt%shape, & + idt%longname, input_attr, this%deflate, this%shuffle, & + this%chunk_face, export_pkg%iper, idxaux, this%nc_fname) + + deallocate (dbl3d) + end subroutine export_layer_3d + !> @brief netcdf export an input array !< subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) - use InputOutputModule, only: lowcase - use MemoryHelperModule, only: memPathSeparator class(Mesh2dDisExportType), intent(inout) :: this character(len=*), intent(in) :: pkgtype character(len=*), intent(in) :: pkgname @@ -188,70 +362,60 @@ subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) real(DP), dimension(:), pointer, contiguous :: dbl1d real(DP), dimension(:, :), pointer, contiguous :: dbl2d real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d - character(len=LINELENGTH) :: pname, vname, nc_varname, gridmap, input_attr - ! - ! -- set package name - pname = pkgname - call lowcase(pname) - ! -- set variable name - vname = idt%tagname - call lowcase(vname) - ! -- set variable name written to file - nc_varname = trim(pname)//'_'//trim(vname) + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: iper, iaux ! - if (this%ogc_wkt /= '') then - ! -- set gridmap variable name - gridmap = this%gridmap_name - else - gridmap = '' - end if + iper = 0 + iaux = 0 ! - if (this%input_attr > 0) then - input_attr = trim(this%modelname)//memPathSeparator//trim(pkgname)// & - memPathSeparator//trim(idt%mf6varname) - else - input_attr = '' - end if + ! -- set package base name + nc_varname = trim(pkgname)//'_'//trim(idt%mf6varname) + ! -- put input attributes + input_attr = this%input_attribute(pkgname, idt) ! select case (idt%datatype) case ('INTEGER1D') call mem_setptr(int1d, idt%mf6varname, mempath) call nc_export_int1d(this%ncid, this%dim_ids, this%x_dim, this%y_dim, & this%var_ids, this%dis, int1d, nc_varname, pkgname, & - idt%tagname, gridmap, idt%shape, idt%longname, & - input_attr, this%deflate, this%shuffle, & - this%chunk_face, this%nc_fname) + idt%tagname, this%gridmap_name, idt%shape, & + idt%longname, input_attr, this%deflate, this%shuffle, & + this%chunk_face, iper, this%nc_fname) case ('INTEGER2D') call mem_setptr(int2d, idt%mf6varname, mempath) call nc_export_int2d(this%ncid, this%dim_ids, this%var_ids, this%dis, & - int2d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + int2d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, this%nc_fname) case ('INTEGER3D') call mem_setptr(int3d, idt%mf6varname, mempath) call nc_export_int3d(this%ncid, this%dim_ids, this%var_ids, this%dis, & - int3d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + int3d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, this%nc_fname) case ('DOUBLE1D') call mem_setptr(dbl1d, idt%mf6varname, mempath) call nc_export_dbl1d(this%ncid, this%dim_ids, this%x_dim, this%y_dim, & this%var_ids, this%dis, dbl1d, nc_varname, pkgname, & - idt%tagname, gridmap, idt%shape, idt%longname, & - input_attr, this%deflate, this%shuffle, & + idt%tagname, this%gridmap_name, idt%shape, & + idt%longname, input_attr, this%deflate, this%shuffle, & this%chunk_face, this%nc_fname) case ('DOUBLE2D') call mem_setptr(dbl2d, idt%mf6varname, mempath) call nc_export_dbl2d(this%ncid, this%dim_ids, this%var_ids, this%dis, & - dbl2d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + dbl2d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, this%nc_fname) case ('DOUBLE3D') call mem_setptr(dbl3d, idt%mf6varname, mempath) call nc_export_dbl3d(this%ncid, this%dim_ids, this%var_ids, this%dis, & - dbl3d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + dbl3d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, iper, iaux, this%nc_fname) case default ! -- no-op, no other datatypes exported end select @@ -437,8 +601,7 @@ end subroutine add_mesh_data subroutine nc_export_int1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & nc_varname, pkgname, tagname, gridmap_name, & shapestr, longname, nc_tag, deflate, shuffle, & - chunk_face, nc_fname) - use InputOutputModule, only: lowcase + chunk_face, iper, nc_fname) integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids integer(I4B), intent(in) :: x_dim @@ -456,6 +619,7 @@ subroutine nc_export_int1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & integer(I4B), intent(in) :: deflate integer(I4B), intent(in) :: shuffle integer(I4B), intent(in) :: chunk_face + integer(I4B), intent(in) :: iper character(len=*), intent(in) :: nc_fname ! -- local integer(I4B), dimension(3) :: dis_shape @@ -467,44 +631,41 @@ subroutine nc_export_int1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & character(len=LINELENGTH) :: longname_l, varname_l ! if (shapestr == 'NROW' .or. & - shapestr == 'NCOL') then + shapestr == 'NCOL' .or. & + shapestr == 'NCPL') then ! - if (shapestr == 'NROW') then + select case (shapestr) + case ('NROW') axis_dim = y_dim - else + case ('NCOL') axis_dim = x_dim - end if - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + case ('NCPL') + axis_dim = dim_ids%nmesh_face + end select + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=0, iper=iper) + longname_l = export_longname(longname, pkgname, tagname, layer=0, iper=iper) ! allocate (var_id(1)) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_INT, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/axis_dim/), var_id(1)), & nc_fname) ! ! -- NROW/NCOL shapes use default chunking - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(1), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(1), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(1), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(1), 'long_name', & longname_l), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add mf6 attr + call ncvar_mf6attr(ncid, var_id(1), 0, iper, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -517,46 +678,30 @@ subroutine nc_export_int1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k, iper=iper) + longname_l = export_longname(longname, pkgname, tagname, layer=k, & + iper=iper) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if - ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) + ! -- defalte and shuffle + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, iper, 0, nc_tag, nc_fname) end do ! ! -- reshape input @@ -604,48 +749,32 @@ subroutine nc_export_int2d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B) :: var_id integer(I4B), dimension(:), pointer, contiguous :: int1d integer(I4B), dimension(1) :: layer_shape - character(len=LINELENGTH) :: longname_l + character(len=LINELENGTH) :: longname_l, varname_l ! - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + ! -- set names + varname_l = export_varname(nc_varname) + longname_l = export_longname(longname, pkgname, tagname, 0) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_INT, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/dim_ids%nmesh_face/), var_id), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id, chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -659,7 +788,6 @@ end subroutine nc_export_int2d subroutine nc_export_int3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, gridmap_name, shapestr, longname, & nc_tag, deflate, shuffle, chunk_face, nc_fname) - use InputOutputModule, only: lowcase integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids type(MeshNCVarIdType), intent(inout) :: var_ids @@ -688,46 +816,29 @@ subroutine nc_export_int3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k) + longname_l = export_longname(longname, pkgname, tagname, k) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, 0, 0, nc_tag, nc_fname) end do ! ! -- exit define mode and write data @@ -748,7 +859,6 @@ subroutine nc_export_dbl1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & nc_varname, pkgname, tagname, gridmap_name, & shapestr, longname, nc_tag, deflate, shuffle, & chunk_face, nc_fname) - use InputOutputModule, only: lowcase integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids integer(I4B), intent(in) :: x_dim @@ -777,44 +887,41 @@ subroutine nc_export_dbl1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & character(len=LINELENGTH) :: longname_l, varname_l ! if (shapestr == 'NROW' .or. & - shapestr == 'NCOL') then + shapestr == 'NCOL') then ! .or. & + !shapestr == 'NCPL') then ! - if (shapestr == 'NROW') then + select case (shapestr) + case ('NROW') axis_dim = y_dim - else + case ('NCOL') axis_dim = x_dim - end if - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + !case ('NCPL') + ! axis_dim = dim_ids%nmesh_face + end select + ! + ! -- set names + varname_l = export_varname(nc_varname) + longname_l = export_longname(longname, pkgname, tagname, 0) ! allocate (var_id(1)) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_DOUBLE, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/axis_dim/), var_id(1)), & nc_fname) ! ! -- NROW/NCOL shapes use default chunking - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(1), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(1), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(1), '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(1), 'long_name', & longname_l), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add mf6 attr + call ncvar_mf6attr(ncid, var_id(1), 0, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -827,46 +934,29 @@ subroutine nc_export_dbl1d(ncid, dim_ids, x_dim, y_dim, var_ids, dis, p_mem, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k) + longname_l = export_longname(longname, pkgname, tagname, k) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if - ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) + ! -- defalte and shuffle + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, 0, 0, nc_tag, nc_fname) end do ! ! -- reshape input @@ -912,50 +1002,34 @@ subroutine nc_export_dbl2d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & character(len=*), intent(in) :: nc_fname ! -- local integer(I4B) :: var_id - character(len=LINELENGTH) :: longname_l + character(len=LINELENGTH) :: longname_l, varname_l real(DP), dimension(:), pointer, contiguous :: dbl1d integer(I4B), dimension(1) :: layer_shape ! - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + ! -- set names + varname_l = export_varname(nc_varname) + longname_l = export_longname(longname, pkgname, tagname, 0) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_DOUBLE, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/dim_ids%nmesh_face/), var_id), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id, chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -968,8 +1042,9 @@ end subroutine nc_export_dbl2d !< subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, gridmap_name, shapestr, longname, & - nc_tag, deflate, shuffle, chunk_face, nc_fname) - use InputOutputModule, only: lowcase + nc_tag, deflate, shuffle, chunk_face, iper, iaux, & + nc_fname) + use ConstantsModule, only: DNODATA integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids type(MeshNCVarIdType), intent(inout) :: var_ids @@ -985,6 +1060,8 @@ subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B), intent(in) :: deflate integer(I4B), intent(in) :: shuffle integer(I4B), intent(in) :: chunk_face + integer(I4B), intent(in) :: iper + integer(I4B), intent(in) :: iaux character(len=*), intent(in) :: nc_fname ! -- local integer(I4B), dimension(:), allocatable :: var_id @@ -992,52 +1069,43 @@ subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & character(len=LINELENGTH) :: longname_l, varname_l integer(I4B), dimension(1) :: layer_shape integer(I4B) :: k + real(DP) :: fill_value + ! + if (iper > 0) then + fill_value = DNODATA + else + fill_value = NF90_FILL_DOUBLE + end if ! allocate (var_id(dis%nlay)) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k, iper=iper, iaux=iaux) + longname_l = export_longname(longname, pkgname, tagname, layer=k, iper=iper) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & - (/NF90_FILL_DOUBLE/)), nc_fname) + (/fill_value/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, iper, iaux, nc_tag, nc_fname) + !end if end do ! ! -- exit define mode and write data diff --git a/src/Utilities/Export/DisNCStructured.f90 b/src/Utilities/Export/DisNCStructured.f90 index 7c52161eb45..328f38b4130 100644 --- a/src/Utilities/Export/DisNCStructured.f90 +++ b/src/Utilities/Export/DisNCStructured.f90 @@ -9,13 +9,13 @@ module DisNCStructuredModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENBIGLINE, LENCOMPONENTNAME, & - LENMEMPATH + LENMEMPATH, LENVARNAME use SimVariablesModule, only: errmsg use SimModule, only: store_error, store_error_filename use MemoryManagerModule, only: mem_setptr use InputDefinitionModule, only: InputParamDefinitionType use CharacterStringModule, only: CharacterStringType - use NCModelExportModule, only: NCBaseModelExportType + use NCModelExportModule, only: NCBaseModelExportType, export_longname use DisModule, only: DisType use NetCDFCommonModule, only: nf_verify use netcdf @@ -28,6 +28,7 @@ module DisNCStructuredModule integer(I4B) :: x !< number of columns integer(I4B) :: y !< number of rows integer(I4B) :: z !< number of layers + integer(I4B) :: ncpl !< number of cells in layer integer(I4B) :: time !< number of steps integer(I4B) :: bnd !< number in boundary contains @@ -66,6 +67,9 @@ module DisNCStructuredModule procedure :: step procedure :: export_input_array procedure :: export_input_arrays + procedure :: package_step_ilayer + procedure :: package_step + procedure :: export_layer_3d procedure :: add_pkg_data procedure :: add_global_att procedure :: define_dim @@ -249,8 +253,6 @@ end subroutine step !> @brief netcdf export an input array !< subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) - use InputOutputModule, only: lowcase - use MemoryHelperModule, only: memPathSeparator class(DisNCStructuredType), intent(inout) :: this character(len=*), intent(in) :: pkgtype character(len=*), intent(in) :: pkgname @@ -262,76 +264,66 @@ subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) real(DP), dimension(:), pointer, contiguous :: dbl1d real(DP), dimension(:, :), pointer, contiguous :: dbl2d real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d - character(len=LINELENGTH) :: pname, vname, nc_varname, gridmap, input_attr + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: iper, iaux ! - ! -- set variable name written to file - pname = pkgname - vname = idt%tagname - call lowcase(pname) - call lowcase(vname) - nc_varname = trim(pname)//'_'//trim(vname) - if (this%ogc_wkt /= '') then - gridmap = this%gridmap_name - else - gridmap = '' - end if + ! -- initialize + iper = 0 + iaux = 0 ! - if (this%input_attr > 0) then - input_attr = trim(this%modelname)//memPathSeparator//trim(pkgname)// & - memPathSeparator//trim(idt%mf6varname) - else - input_attr = '' - end if + ! -- set variable name and input attribute string + nc_varname = export_varname(pkgname, idt) + input_attr = this%input_attribute(pkgname, idt) ! select case (idt%datatype) case ('INTEGER1D') call mem_setptr(int1d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & int1d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, iper, this%nc_fname) case ('INTEGER2D') call mem_setptr(int2d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & int2d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, this%nc_fname) case ('INTEGER3D') call mem_setptr(int3d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & int3d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, this%nc_fname) case ('DOUBLE1D') call mem_setptr(dbl1d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & dbl1d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, iper, this%nc_fname) case ('DOUBLE2D') call mem_setptr(dbl2d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & dbl2d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, this%nc_fname) case ('DOUBLE3D') call mem_setptr(dbl3d, idt%mf6varname, mempath) call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & dbl3d, nc_varname, pkgname, idt%tagname, & - idt%shape, idt%longname, input_attr, gridmap, & - this%latlon, this%deflate, this%shuffle, & - this%chunk_z, this%chunk_y, this%chunk_x, & - this%nc_fname) + idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, iper, iaux, this%nc_fname) case default ! -- no-op, no other datatypes exported end select @@ -365,6 +357,232 @@ subroutine export_input_arrays(this, pkgtype, pkgname, mempath, param_dfns) end do end subroutine export_input_arrays + !> @brief netcdf export package dynamic input with ilayer index variable + !< + subroutine package_step_ilayer(this, export_pkg, ilayer_varname, ilayer) + use ConstantsModule, only: DNODATA, DZERO + use TdisModule, only: kper + use NCModelExportModule, only: ExportPackageType + use DefinitionSelectModule, only: get_param_definition_type + class(DisNCStructuredType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + character(len=*), intent(in) :: ilayer_varname + integer(I4B), intent(in) :: ilayer + type(InputParamDefinitionType), pointer :: idt + ! -- local + integer(I4B), dimension(:), pointer, contiguous :: int1d + real(DP), dimension(:), pointer, contiguous :: dbl1d + real(DP), dimension(:, :), pointer, contiguous :: dbl2d + integer(I4B), dimension(:), pointer, contiguous :: ialayer + real(DP), dimension(:), contiguous, pointer :: dbl1d_ptr + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: n, iparam, nvals + logical(LGP) :: ilayer_read + ! + ! -- initialize + nullify (ialayer) + ilayer_read = .false. + ! + ! -- set pointer to ilayer variable + call mem_setptr(ialayer, export_pkg%param_names(ilayer), & + export_pkg%mf6_input%mempath) + ! + ! -- check if layer index variable was read + if (export_pkg%param_reads(ilayer)%invar == 1) then + ilayer_read = .true. + end if + ! + ! -- export defined period input + do iparam = 1, export_pkg%nparam + ! + ! -- check if variable was read this period + if (export_pkg%param_reads(iparam)%invar < 1) cycle + ! + ! -- set input definition + idt => & + get_param_definition_type(export_pkg%mf6_input%param_dfns, & + export_pkg%mf6_input%component_type, & + export_pkg%mf6_input%subcomponent_type, & + 'PERIOD', export_pkg%param_names(iparam), & + this%nc_fname) + ! + ! -- set variable name and input attrs + nc_varname = export_varname(export_pkg%mf6_input%subcomponent_name, idt, & + iper=kper) + input_attr = this%input_attribute(export_pkg%mf6_input%subcomponent_name, & + idt) + ! + ! -- export arrays + select case (idt%datatype) + case ('INTEGER1D') + call mem_setptr(int1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & + int1d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, & + idt%tagname, idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, export_pkg%iper, this%nc_fname) + case ('DOUBLE1D') + call mem_setptr(dbl1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call this%export_layer_3d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr) + case ('DOUBLE2D') + call mem_setptr(dbl2d, idt%mf6varname, export_pkg%mf6_input%mempath) + nvals = this%dis%ncol * this%dis%nrow + ! + do n = 1, size(dbl2d, dim=1) ! naux + dbl1d_ptr(1:nvals) => dbl2d(n, :) + if (all(dbl1d_ptr == DZERO)) then + else + call this%export_layer_3d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d_ptr, nc_varname, input_attr, n) + end if + end do + case default + ! + errmsg = 'EXPORT ilayer unsupported datatype='//trim(idt%datatype) + call store_error(errmsg, .true.) + end select + end do + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + ! + ! -- return + return + end subroutine package_step_ilayer + + !> @brief netcdf export package dynamic input + !< + subroutine package_step(this, export_pkg) + use TdisModule, only: kper + use NCModelExportModule, only: ExportPackageType + use DefinitionSelectModule, only: get_param_definition_type + class(DisNCStructuredType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + integer(I4B), dimension(:), pointer, contiguous :: int1d + real(DP), dimension(:), pointer, contiguous :: dbl1d + type(InputParamDefinitionType), pointer :: idt + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: iparam + ! + do iparam = 1, export_pkg%nparam + ! + ! -- set input definition + idt => get_param_definition_type(export_pkg%mf6_input%param_dfns, & + export_pkg%mf6_input%component_type, & + export_pkg%mf6_input%subcomponent_type, & + 'PERIOD', export_pkg%param_names(iparam), & + this%nc_fname) + ! + ! -- set variable name and input attribute string + nc_varname = export_varname(export_pkg%mf6_input%subcomponent_name, idt, & + iper=kper) + input_attr = this%input_attribute(export_pkg%mf6_input%subcomponent_name, & + idt) + ! + ! -- export arrays + select case (idt%datatype) + case ('INTEGER1D') + ! + call mem_setptr(int1d, export_pkg%param_names(iparam), & + export_pkg%mf6_input%mempath) + call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & + int1d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, & + idt%tagname, idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, kper, this%nc_fname) + ! + case ('DOUBLE1D') + ! + call mem_setptr(dbl1d, export_pkg%param_names(iparam), & + export_pkg%mf6_input%mempath) + call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, & + dbl1d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, & + idt%tagname, idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, & + this%chunk_x, kper, this%nc_fname) + ! + case default + ! + errmsg = 'EXPORT unsupported datatype='//trim(idt%datatype) + call store_error(errmsg, .true.) + ! + end select + end do + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + end subroutine package_step + + !> @brief export layer variable as full grid + !< + subroutine export_layer_3d(this, export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr, iaux) + use ConstantsModule, only: DNODATA, DZERO + use TdisModule, only: kper + use NCModelExportModule, only: ExportPackageType + class(DisNCStructuredType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + type(InputParamDefinitionType), pointer, intent(in) :: idt + logical(LGP), intent(in) :: ilayer_read + integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ialayer + real(DP), dimension(:), pointer, contiguous, intent(in) :: dbl1d + character(len=*), intent(inout) :: nc_varname + character(len=*), intent(in) :: input_attr + integer(I4B), optional, intent(in) :: iaux + ! -- local + real(DP), dimension(:, :, :), pointer, contiguous :: dbl3d + integer(I4B) :: n, i, j, k, nvals, idxaux + real(DP), dimension(:, :), contiguous, pointer :: dbl2d_ptr + ! + ! -- initialize + idxaux = 0 + if (present(iaux)) then + nc_varname = export_varname(export_pkg%mf6_input%subcomponent_name, & + idt, iper=kper, iaux=iaux) + idxaux = iaux + end if + + allocate (dbl3d(export_pkg%mshape(3), export_pkg%mshape(2), & + export_pkg%mshape(1))) + ! + if (ilayer_read) then + do k = 1, size(dbl3d, dim=3) + n = 0 + do i = 1, size(dbl3d, dim=2) + do j = 1, size(dbl3d, dim=1) + n = n + 1 + if (ialayer(n) == k) then + dbl3d(j, i, k) = dbl1d(n) + else + dbl3d(j, i, k) = DNODATA + end if + end do + end do + end do + else + dbl3d = DNODATA + nvals = export_pkg%mshape(3) * export_pkg%mshape(2) + dbl2d_ptr(1:export_pkg%mshape(3), 1:export_pkg%mshape(2)) => dbl1d(1:nvals) + dbl3d(:, :, 1) = dbl2d_ptr(:, :) + end if + ! + call nc_export_array(this%ncid, this%dim_ids, this%var_ids, this%dis, dbl3d, & + nc_varname, export_pkg%mf6_input%subcomponent_name, & + idt%tagname, idt%shape, idt%longname, input_attr, & + this%gridmap_name, this%latlon, this%deflate, & + this%shuffle, this%chunk_z, this%chunk_y, this%chunk_x, & + export_pkg%iper, idxaux, this%nc_fname) + + deallocate (dbl3d) + end subroutine export_layer_3d + !> @brief determine packages to write gridded input !< subroutine add_pkg_data(this) @@ -540,6 +758,12 @@ subroutine define_dim(this) call nf_verify(nf90_def_var(this%ncid, 'x_bnds', NF90_DOUBLE, & (/this%dim_ids%bnd, this%dim_ids%x/), & this%var_ids%x_bnds), this%nc_fname) + + ! + ! -- NCPL dimension + call nf_verify(nf90_def_dim(this%ncid, 'ncpl', & + this%dis%ncol * this%dis%nrow, & + this%dim_ids%ncpl), this%nc_fname) end subroutine define_dim !> @brief create the model layer dependent variables @@ -563,14 +787,10 @@ subroutine define_dependent(this) this%nc_fname) end if ! -- deflate and shuffle - if (this%deflate >= 0) then - call nf_verify(nf90_def_var_deflate(this%ncid, & - this%var_ids%dependent, & - shuffle=this%shuffle, deflate=1, & - deflate_level=this%deflate), & - this%nc_fname) - end if + call ncvar_deflate(this%ncid, this%var_ids%dependent, this%deflate, & + this%shuffle, this%nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, & 'units', 'm'), this%nc_fname) call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, & @@ -580,19 +800,10 @@ subroutine define_dependent(this) this%annotation%longname), this%nc_fname) call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, '_FillValue', & (/DHNOFLO/)), this%nc_fname) - if (this%ogc_wkt /= '') then - if (this%latlon) then - call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, & - 'coordinates', 'lon lat'), & - this%nc_fname) - else - call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, & - 'coordinates', 'x y'), this%nc_fname) - end if - call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent, & - 'grid_mapping', this%gridmap_name), & - this%nc_fname) - end if + ! + ! -- add grid mapping + call ncvar_gridmap(this%ncid, this%var_ids%dependent, this%gridmap_name, & + this%latlon, this%nc_fname) end subroutine define_dependent !> @brief create the file grid mapping container variable @@ -728,13 +939,107 @@ subroutine add_grid_data(this) deallocate (dbl2d) end subroutine add_grid_data + !> @brief define 2d variable chunking + !< + subroutine ncvar_chunk2d(ncid, varid, chunk_x, chunk_y, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: chunk_x + integer(I4B), intent(in) :: chunk_y + character(len=*), intent(in) :: nc_fname + ! + if (chunk_y > 0 .and. chunk_x > 0) then + call nf_verify(nf90_def_var_chunking(ncid, varid, NF90_CHUNKED, & + (/chunk_x, chunk_y/)), nc_fname) + end if + end subroutine ncvar_chunk2d + + !> @brief define 3d variable chunking + !< + subroutine ncvar_chunk3d(ncid, varid, chunk_x, chunk_y, chunk_z, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: chunk_x + integer(I4B), intent(in) :: chunk_y + integer(I4B), intent(in) :: chunk_z + character(len=*), intent(in) :: nc_fname + ! + if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then + call nf_verify(nf90_def_var_chunking(ncid, varid, NF90_CHUNKED, & + (/chunk_x, chunk_y, chunk_z/)), & + nc_fname) + end if + end subroutine ncvar_chunk3d + + !> @brief define variable compression + !< + subroutine ncvar_deflate(ncid, varid, deflate, shuffle, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: deflate + integer(I4B), intent(in) :: shuffle + character(len=*), intent(in) :: nc_fname + ! -- deflate and shuffle + if (deflate >= 0) then + call nf_verify(nf90_def_var_deflate(ncid, varid, shuffle=shuffle, & + deflate=1, deflate_level=deflate), & + nc_fname) + end if + end subroutine ncvar_deflate + + !> @brief put variable gridmap attributes + !< + subroutine ncvar_gridmap(ncid, varid, gridmap_name, latlon, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + character(len=*), intent(in) :: gridmap_name + logical(LGP), intent(in) :: latlon + character(len=*), intent(in) :: nc_fname + ! + if (gridmap_name /= '') then + if (latlon) then + call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'lon lat'), & + nc_fname) + else + call nf_verify(nf90_put_att(ncid, varid, 'coordinates', 'x y'), & + nc_fname) + end if + call nf_verify(nf90_put_att(ncid, varid, 'grid_mapping', gridmap_name), & + nc_fname) + end if + end subroutine ncvar_gridmap + + !> @brief put variable internal modflow6 attributes + !< + subroutine ncvar_mf6attr(ncid, varid, iper, iaux, nc_tag, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: iper + integer(I4B), intent(in) :: iaux + character(len=*), intent(in) :: nc_tag + character(len=*), intent(in) :: nc_fname + ! + if (nc_tag /= '') then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_input', & + nc_tag), nc_fname) + if (iper > 0) then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_iper', & + iper), nc_fname) + end if + ! + if (iaux > 0) then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_iaux', & + iaux), nc_fname) + end if + end if + end subroutine ncvar_mf6attr + !> @brief netcdf export 1D integer !< subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, shapestr, longname, nc_tag, & gridmap_name, latlon, deflate, shuffle, chunk_z, & - chunk_y, chunk_x, nc_fname) - use InputOutputModule, only: lowcase + chunk_y, chunk_x, iper, nc_fname) integer(I4B), intent(in) :: ncid type(StructuredNCDimIdType), intent(inout) :: dim_ids type(StructuredNCVarIdType), intent(inout) :: var_ids @@ -753,17 +1058,26 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B), intent(in) :: chunk_z integer(I4B), intent(in) :: chunk_y integer(I4B), intent(in) :: chunk_x + integer(I4B), intent(in) :: iper character(len=*), intent(in) :: nc_fname ! -- local integer(I4B) :: var_id, axis_sz + character(len=LINELENGTH) :: longname_l ! if (shapestr == 'NROW' .or. & - shapestr == 'NCOL') then - if (shapestr == 'NROW') then + shapestr == 'NCOL' .or. & + shapestr == 'NCPL') then + ! + select case (shapestr) + case ('NROW') axis_sz = dim_ids%y - else + case ('NCOL') axis_sz = dim_ids%x - end if + case ('NCPL') + axis_sz = dim_ids%ncpl + end select + ! + longname_l = export_longname(longname, pkgname, tagname, layer=0, iper=iper) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) @@ -772,21 +1086,16 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- NROW/NCOL shapes use default chunking - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & - longname), nc_fname) + longname_l), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add mf6 attr + call ncvar_mf6attr(ncid, var_id, iper, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -802,38 +1111,19 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y, chunk_z/)), & - nc_fname) - end if + call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -878,37 +1168,19 @@ subroutine nc_export_int2d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y/)), nc_fname) - end if + call ncvar_chunk2d(ncid, var_id, chunk_x, chunk_y, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -923,7 +1195,6 @@ subroutine nc_export_int3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, shapestr, longname, nc_tag, & gridmap_name, latlon, deflate, shuffle, chunk_z, & chunk_y, chunk_x, nc_fname) - use InputOutputModule, only: lowcase integer(I4B), intent(in) :: ncid type(StructuredNCDimIdType), intent(inout) :: dim_ids type(StructuredNCVarIdType), intent(inout) :: var_ids @@ -953,38 +1224,19 @@ subroutine nc_export_int3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y, chunk_z/)), & - nc_fname) - end if + call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -998,8 +1250,8 @@ end subroutine nc_export_int3d subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, shapestr, longname, nc_tag, & gridmap_name, latlon, deflate, shuffle, chunk_z, & - chunk_y, chunk_x, nc_fname) - use InputOutputModule, only: lowcase + chunk_y, chunk_x, iper, nc_fname) + use ConstantsModule, only: DNODATA integer(I4B), intent(in) :: ncid type(StructuredNCDimIdType), intent(inout) :: dim_ids type(StructuredNCVarIdType), intent(inout) :: var_ids @@ -1018,17 +1270,25 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B), intent(in) :: chunk_z integer(I4B), intent(in) :: chunk_y integer(I4B), intent(in) :: chunk_x + integer(I4B), intent(in) :: iper character(len=*), intent(in) :: nc_fname ! -- local integer(I4B) :: var_id, axis_sz + real(DP) :: fill_value + character(len=LINELENGTH) :: longname_l ! if (shapestr == 'NROW' .or. & - shapestr == 'NCOL') then - if (shapestr == 'NROW') then + shapestr == 'NCOL' .or. & + shapestr == 'NCPL') then + ! + select case (shapestr) + case ('NROW') axis_sz = dim_ids%y - else + case ('NCOL') axis_sz = dim_ids%x - end if + case ('NCPL') + axis_sz = dim_ids%ncpl + end select ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) @@ -1037,21 +1297,16 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- NROW/NCOL shapes use default chunking - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add mf6 attr + call ncvar_mf6attr(ncid, var_id, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -1059,6 +1314,14 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) else + ! + if (iper > 0) then + fill_value = DNODATA + else + fill_value = NF90_FILL_DOUBLE + end if + ! + longname_l = export_longname(longname, pkgname, tagname, layer=0, iper=iper) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) @@ -1067,38 +1330,19 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y, chunk_z/)), & - nc_fname) - end if + call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & - (/NF90_FILL_DOUBLE/)), nc_fname) + (/fill_value/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & - longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if + longname_l), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, iper, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -1143,37 +1387,19 @@ subroutine nc_export_dbl2d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y/)), nc_fname) - end if + call ncvar_chunk2d(ncid, var_id, chunk_x, chunk_y, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -1187,8 +1413,8 @@ end subroutine nc_export_dbl2d subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, shapestr, longname, nc_tag, & gridmap_name, latlon, deflate, shuffle, chunk_z, & - chunk_y, chunk_x, nc_fname) - use InputOutputModule, only: lowcase + chunk_y, chunk_x, iper, iaux, nc_fname) + use ConstantsModule, only: DNODATA integer(I4B), intent(in) :: ncid type(StructuredNCDimIdType), intent(inout) :: dim_ids type(StructuredNCVarIdType), intent(inout) :: var_ids @@ -1207,9 +1433,21 @@ subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B), intent(in) :: chunk_z integer(I4B), intent(in) :: chunk_y integer(I4B), intent(in) :: chunk_x + integer(I4B), intent(in) :: iper + integer(I4B), intent(in) :: iaux character(len=*), intent(in) :: nc_fname ! -- local integer(I4B) :: var_id + real(DP) :: fill_value + character(len=LINELENGTH) :: longname_l + ! + if (iper > 0) then + fill_value = DNODATA + else + fill_value = NF90_FILL_DOUBLE + end if + ! + longname_l = export_longname(longname, pkgname, tagname, layer=0, iper=iper) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) @@ -1218,38 +1456,19 @@ subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) ! ! -- apply chunking parameters - if (chunk_z > 0 .and. chunk_y > 0 .and. chunk_x > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id, NF90_CHUNKED, & - (/chunk_x, chunk_y, chunk_z/)), & - nc_fname) - end if + call ncvar_chunk3d(ncid, var_id, chunk_x, chunk_y, chunk_z, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id, shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id, deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id, '_FillValue', & - (/NF90_FILL_DOUBLE/)), nc_fname) + (/fill_value/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id, 'long_name', & - longname), nc_fname) - if (gridmap_name /= '') then - if (latlon) then - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'lon lat'), & - nc_fname) - else - call nf_verify(nf90_put_att(ncid, var_id, 'coordinates', 'x y'), & - nc_fname) - end if - call nf_verify(nf90_put_att(ncid, var_id, 'grid_mapping', gridmap_name), & - nc_fname) - end if + longname_l), nc_fname) ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id, 'modflow6_input', & - nc_tag), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id, gridmap_name, latlon, nc_fname) + call ncvar_mf6attr(ncid, var_id, iper, iaux, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -1258,4 +1477,31 @@ subroutine nc_export_dbl3d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & nc_fname) end subroutine nc_export_dbl3d + !> @brief build netcdf variable name + !< + function export_varname(pkgname, idt, iper, iaux) result(varname) + use InputOutputModule, only: lowcase + character(len=*), intent(in) :: pkgname + type(InputParamDefinitionType), pointer, intent(in) :: idt + integer(I4B), optional, intent(in) :: iper + integer(I4B), optional, intent(in) :: iaux + character(len=LINELENGTH) :: varname + character(len=LINELENGTH) :: pname, vname + ! + pname = pkgname + vname = idt%mf6varname + call lowcase(pname) + call lowcase(vname) + if (present(iper)) then + if (present(iaux)) then + write (varname, '(a,i0,a,i0)') trim(pname)//'_'//trim(vname)// & + '_p', iper, 'a', iaux + else + write (varname, '(a,i0)') trim(pname)//'_'//trim(vname)//'_p', iper + end if + else + varname = trim(pname)//'_'//trim(vname) + end if + end function export_varname + end module DisNCStructuredModule diff --git a/src/Utilities/Export/DisvNCMesh.f90 b/src/Utilities/Export/DisvNCMesh.f90 index 87c6f024e18..92655c7b380 100644 --- a/src/Utilities/Export/DisvNCMesh.f90 +++ b/src/Utilities/Export/DisvNCMesh.f90 @@ -10,10 +10,15 @@ module MeshDisvModelModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENBIGLINE, LENCOMPONENTNAME, & LENMEMPATH + use SimVariablesModule, only: errmsg + use SimModule, only: store_error, store_error_filename use MemoryManagerModule, only: mem_setptr use InputDefinitionModule, only: InputParamDefinitionType use CharacterStringModule, only: CharacterStringType - use MeshModelModule, only: Mesh2dModelType, MeshNCDimIdType, MeshNCVarIdType + use MeshModelModule, only: Mesh2dModelType, MeshNCDimIdType, MeshNCVarIdType, & + ncvar_chunk, ncvar_deflate, ncvar_gridmap, & + ncvar_mf6attr, export_varname + use NCModelExportModule, only: export_longname use DisvModule, only: DisvType use NetCDFCommonModule, only: nf_verify use netcdf @@ -31,6 +36,9 @@ module MeshDisvModelModule procedure :: df procedure :: step procedure :: export_input_array + procedure :: package_step_ilayer + procedure :: package_step + procedure :: export_layer_2d procedure :: define_dim procedure :: add_mesh_data end type Mesh2dDisvExportType @@ -170,11 +178,172 @@ subroutine step(this) nullify (dbl2d) end subroutine step + !> @brief netcdf export package dynamic input with ilayer index variable + !< + subroutine package_step_ilayer(this, export_pkg, ilayer_varname, ilayer) + use ConstantsModule, only: DNODATA, DZERO + use TdisModule, only: kper + use DefinitionSelectModule, only: get_param_definition_type + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisvExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + character(len=*), intent(in) :: ilayer_varname + integer(I4B), intent(in) :: ilayer + ! -- local + type(InputParamDefinitionType), pointer :: idt + integer(I4B), dimension(:), pointer, contiguous :: int1d + real(DP), dimension(:), pointer, contiguous :: dbl1d + real(DP), dimension(:, :), pointer, contiguous :: dbl2d + integer(I4B), dimension(:), pointer, contiguous :: ialayer + real(DP), dimension(:), contiguous, pointer :: dbl1d_ptr + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: n, iparam, nvals + logical(LGP) :: ilayer_read + ! + ! -- initialize + nullify (ialayer) + ilayer_read = .false. + ! + ! -- set pointer to ilayer variable + call mem_setptr(ialayer, export_pkg%param_names(ilayer), & + export_pkg%mf6_input%mempath) + ! + ! -- check if layer index variable was read + if (export_pkg%param_reads(ilayer)%invar == 1) then + ilayer_read = .true. + end if + ! + ! -- export defined period input + do iparam = 1, export_pkg%nparam + ! + ! -- check if variable was read this period + if (export_pkg%param_reads(iparam)%invar < 1) cycle + ! + ! -- set input definition + idt => & + get_param_definition_type(export_pkg%mf6_input%param_dfns, & + export_pkg%mf6_input%component_type, & + export_pkg%mf6_input%subcomponent_type, & + 'PERIOD', export_pkg%param_names(iparam), '') + ! + ! -- set variable name and input string + nc_varname = trim(export_pkg%mf6_input%subcomponent_name)//'_'// & + trim(idt%mf6varname) + input_attr = this%input_attribute(export_pkg%mf6_input%subcomponent_name, & + idt) + ! + ! -- export arrays + select case (idt%datatype) + case ('INTEGER1D') + call mem_setptr(int1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call nc_export_int1d(this%ncid, this%dim_ids, this%var_ids, this%disv, & + int1d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, & + idt%tagname, this%gridmap_name, idt%shape, & + idt%longname, input_attr, this%deflate, & + this%shuffle, this%chunk_face, kper, this%nc_fname) + case ('DOUBLE1D') + call mem_setptr(dbl1d, idt%mf6varname, export_pkg%mf6_input%mempath) + call this%export_layer_2d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr) + case ('DOUBLE2D') + call mem_setptr(dbl2d, idt%mf6varname, export_pkg%mf6_input%mempath) + nvals = this%disv%ncpl + ! + do n = 1, size(dbl2d, dim=1) !naux + dbl1d_ptr(1:nvals) => dbl2d(n, :) + if (all(dbl1d_ptr == DZERO)) then + else + call this%export_layer_2d(export_pkg, idt, ilayer_read, ialayer, & + dbl1d_ptr, nc_varname, input_attr, n) + end if + end do + case default + ! + errmsg = 'EXPORT ilayer unsupported datatype='//trim(idt%datatype) + call store_error(errmsg, .true.) + end select + end do + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + ! + ! -- return + return + end subroutine package_step_ilayer + + !> @brief netcdf export package dynamic input + !< + subroutine package_step(this, export_pkg) + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisvExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + errmsg = 'NetCDF period export not supported for model='// & + trim(this%modelname)//', package='// & + trim(export_pkg%mf6_input%subcomponent_name) + call store_error(errmsg, .true.) + ! + ! -- synchronize file + call nf_verify(nf90_sync(this%ncid), this%nc_fname) + end subroutine package_step + + !> @brief export layer variable as full grid + !< + subroutine export_layer_2d(this, export_pkg, idt, ilayer_read, ialayer, & + dbl1d, nc_varname, input_attr, iaux) + use ConstantsModule, only: DNODATA, DZERO + use NCModelExportModule, only: ExportPackageType + class(Mesh2dDisvExportType), intent(inout) :: this + class(ExportPackageType), pointer, intent(in) :: export_pkg + type(InputParamDefinitionType), pointer, intent(in) :: idt + logical(LGP), intent(in) :: ilayer_read + integer(I4B), dimension(:), pointer, contiguous, intent(in) :: ialayer + real(DP), dimension(:), pointer, contiguous, intent(in) :: dbl1d + character(len=*), intent(in) :: nc_varname + character(len=*), intent(in) :: input_attr + integer(I4B), optional, intent(in) :: iaux + ! -- local + real(DP), dimension(:, :), pointer, contiguous :: dbl2d + integer(I4B) :: n, j, k, idxaux + ! + ! -- initialize + idxaux = 0 + if (present(iaux)) then + idxaux = iaux + end if + + allocate (dbl2d(export_pkg%mshape(2), export_pkg%mshape(1))) + ! + if (ilayer_read) then + do k = 1, size(dbl2d, dim=2) + n = 0 + do j = 1, size(dbl2d, dim=1) + n = n + 1 + if (ialayer(n) == k) then + dbl2d(j, k) = dbl1d(n) + else + dbl2d(j, k) = DNODATA + end if + end do + end do + else + dbl2d = DNODATA + dbl2d(:, 1) = dbl1d(:) + end if + ! + call nc_export_dbl2d(this%ncid, this%dim_ids, this%var_ids, this%disv, & + dbl2d, nc_varname, & + export_pkg%mf6_input%subcomponent_name, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, input_attr, & + this%deflate, this%shuffle, this%chunk_face, & + export_pkg%iper, idxaux, this%nc_fname) + ! + deallocate (dbl2d) + end subroutine export_layer_2d + !> @brief netcdf export an input array !< subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) - use InputOutputModule, only: lowcase - use MemoryHelperModule, only: memPathSeparator class(Mesh2dDisvExportType), intent(inout) :: this character(len=*), intent(in) :: pkgtype character(len=*), intent(in) :: pkgname @@ -184,56 +353,46 @@ subroutine export_input_array(this, pkgtype, pkgname, mempath, idt) integer(I4B), dimension(:, :), pointer, contiguous :: int2d real(DP), dimension(:), pointer, contiguous :: dbl1d real(DP), dimension(:, :), pointer, contiguous :: dbl2d - character(len=LINELENGTH) :: pname, vname, nc_varname, gridmap, input_attr - ! - ! -- set package name - pname = pkgname - call lowcase(pname) - ! -- set variable name - vname = idt%tagname - call lowcase(vname) - ! -- set variable name written to file - nc_varname = trim(pname)//'_'//trim(vname) - ! - if (this%ogc_wkt /= '') then - ! -- set gridmap variable name - gridmap = this%gridmap_name - else - gridmap = '' - end if + character(len=LINELENGTH) :: nc_varname, input_attr + integer(I4B) :: iper, iaux ! - if (this%input_attr > 0) then - input_attr = trim(this%modelname)//memPathSeparator//trim(pkgname)// & - memPathSeparator//trim(idt%mf6varname) - else - input_attr = '' - end if + iper = 0 + iaux = 0 + ! + ! -- set package base name + nc_varname = trim(pkgname)//'_'//trim(idt%mf6varname) + ! -- put input attributes + input_attr = this%input_attribute(pkgname, idt) ! select case (idt%datatype) case ('INTEGER1D') call mem_setptr(int1d, idt%mf6varname, mempath) call nc_export_int1d(this%ncid, this%dim_ids, this%var_ids, this%disv, & - int1d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + int1d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, iper, this%nc_fname) case ('INTEGER2D') call mem_setptr(int2d, idt%mf6varname, mempath) call nc_export_int2d(this%ncid, this%dim_ids, this%var_ids, this%disv, & - int2d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + int2d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, this%nc_fname) case ('DOUBLE1D') call mem_setptr(dbl1d, idt%mf6varname, mempath) call nc_export_dbl1d(this%ncid, this%dim_ids, this%var_ids, this%disv, & - dbl1d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + dbl1d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, this%nc_fname) case ('DOUBLE2D') call mem_setptr(dbl2d, idt%mf6varname, mempath) call nc_export_dbl2d(this%ncid, this%dim_ids, this%var_ids, this%disv, & - dbl2d, nc_varname, pkgname, idt%tagname, gridmap, & - idt%shape, idt%longname, input_attr, this%deflate, & - this%shuffle, this%chunk_face, this%nc_fname) + dbl2d, nc_varname, pkgname, idt%tagname, & + this%gridmap_name, idt%shape, idt%longname, & + input_attr, this%deflate, this%shuffle, & + this%chunk_face, iper, iaux, this%nc_fname) case default ! -- no-op, no other datatypes exported end select @@ -246,7 +405,6 @@ subroutine define_dim(this) use SimVariablesModule, only: isim_mode class(Mesh2dDisvExportType), intent(inout) :: this integer(I4B), dimension(:), contiguous, pointer :: ncvert - integer(I4B) :: ncpl_dim ! ! -- set pointers to input context call mem_setptr(ncvert, 'NCVERT', this%dis_mempath) @@ -283,8 +441,6 @@ subroutine define_dim(this) ! -- ncpl, nlay call nf_verify(nf90_def_dim(this%ncid, 'nlay', this%disv%nlay, & this%dim_ids%nlay), this%nc_fname) - call nf_verify(nf90_def_dim(this%ncid, 'ncpl', this%disv%ncpl, & - ncpl_dim), this%nc_fname) end subroutine define_dim !> @brief netcdf export add mesh information @@ -390,8 +546,7 @@ end subroutine add_mesh_data !< subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, gridmap_name, shapestr, longname, & - nc_tag, deflate, shuffle, chunk_face, nc_fname) - use InputOutputModule, only: lowcase + nc_tag, deflate, shuffle, chunk_face, iper, nc_fname) integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids type(MeshNCVarIdType), intent(inout) :: var_ids @@ -407,6 +562,7 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & integer(I4B), intent(in) :: deflate integer(I4B), intent(in) :: shuffle integer(I4B), intent(in) :: chunk_face + integer(I4B), intent(in) :: iper character(len=*), intent(in) :: nc_fname ! -- local integer(I4B), dimension(2) :: dis_shape @@ -416,49 +572,34 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & character(len=LINELENGTH) :: longname_l, varname_l ! if (shapestr == 'NCPL') then - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + ! + ! -- set names + varname_l = export_varname(nc_varname) + longname_l = export_longname(longname, pkgname, tagname, layer=0, iper=iper) ! allocate (var_id(1)) axis_sz = dim_ids%nmesh_face ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_INT, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/axis_sz/), var_id(1)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(1), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(1), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(1), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(1), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(1), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(1), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(1), 'grid_mapping', & - gridmap_name), nc_fname) - end if - ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'modflow6_input', & - nc_tag), nc_fname) - end if + ! + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(1), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(1), 0, iper, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -471,46 +612,30 @@ subroutine nc_export_int1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k, iper=iper) + longname_l = export_longname(longname, pkgname, tagname, layer=k, & + iper=iper) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if - ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) + ! -- defalte and shuffle + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, iper, 0, nc_tag, nc_fname) end do ! ! -- reshape input @@ -561,46 +686,29 @@ subroutine nc_export_int2d(ncid, dim_ids, var_ids, disv, p_mem, nc_varname, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, disv%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k) + longname_l = export_longname(longname, pkgname, tagname, k) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_INT, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_INT/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if - ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, 0, 0, nc_tag, nc_fname) end do ! ! -- exit define mode and write data @@ -617,7 +725,6 @@ end subroutine nc_export_int2d subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & pkgname, tagname, gridmap_name, shapestr, longname, & nc_tag, deflate, shuffle, chunk_face, nc_fname) - use InputOutputModule, only: lowcase integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids type(MeshNCVarIdType), intent(inout) :: var_ids @@ -642,49 +749,34 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & character(len=LINELENGTH) :: longname_l, varname_l ! if (shapestr == 'NCPL') then - ! -- set long_name attribute - if (longname /= '') then - longname_l = trim(longname) - else - longname_l = trim(pkgname)//' '//trim(tagname) - end if + ! + ! -- set names + varname_l = export_varname(nc_varname) + longname_l = export_longname(longname, pkgname, tagname, 0) ! allocate (var_id(1)) axis_sz = dim_ids%nmesh_face ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) - call nf_verify(nf90_def_var(ncid, nc_varname, NF90_DOUBLE, & + call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/axis_sz/), var_id(1)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(1), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(1), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(1), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(1), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(1), '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(1), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(1), 'grid_mapping', & - gridmap_name), nc_fname) - end if - ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(1), 'modflow6_input', & - nc_tag), nc_fname) - end if + ! + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(1), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(1), 0, 0, 0, nc_tag, nc_fname) ! ! -- exit define mode and write data call nf_verify(nf90_enddef(ncid), nc_fname) @@ -697,46 +789,29 @@ subroutine nc_export_dbl1d(ncid, dim_ids, var_ids, dis, p_mem, nc_varname, & ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, dis%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0,a)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k) + longname_l = export_longname(longname, pkgname, tagname, k) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & (/NF90_FILL_DOUBLE/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, 0, 0, nc_tag, nc_fname) end do ! ! -- reshape input @@ -760,7 +835,9 @@ end subroutine nc_export_dbl1d !< subroutine nc_export_dbl2d(ncid, dim_ids, var_ids, disv, p_mem, nc_varname, & pkgname, tagname, gridmap_name, shapestr, longname, & - nc_tag, deflate, shuffle, chunk_face, nc_fname) + nc_tag, deflate, shuffle, chunk_face, iper, iaux, & + nc_fname) + use ConstantsModule, only: DNODATA integer(I4B), intent(in) :: ncid type(MeshNCDimIdType), intent(inout) :: dim_ids type(MeshNCVarIdType), intent(inout) :: var_ids @@ -776,57 +853,49 @@ subroutine nc_export_dbl2d(ncid, dim_ids, var_ids, disv, p_mem, nc_varname, & integer(I4B), intent(in) :: deflate integer(I4B), intent(in) :: shuffle integer(I4B), intent(in) :: chunk_face + integer(I4B), intent(in) :: iper + integer(I4B), intent(in) :: iaux character(len=*), intent(in) :: nc_fname ! -- local integer(I4B), dimension(:), allocatable :: var_id character(len=LINELENGTH) :: longname_l, varname_l integer(I4B) :: k + real(DP) :: fill_value + ! + if (iper > 0) then + fill_value = DNODATA + else + fill_value = NF90_FILL_DOUBLE + end if ! allocate (var_id(disv%nlay)) ! ! -- reenter define mode and create variable call nf_verify(nf90_redef(ncid), nc_fname) do k = 1, disv%nlay - write (varname_l, '(a,i0)') trim(nc_varname)//'_l', k - if (longname /= '') then - write (longname_l, '(a,i0,a)') trim(longname)//' (layer ', k, ')' - else - write (longname_l, '(a,i0)') trim(pkgname)//' '//trim(tagname)// & - ' (layer ', k, ')' - end if + ! + ! -- set names + varname_l = export_varname(nc_varname, layer=k, iper=iper, iaux=iaux) + longname_l = export_longname(longname, pkgname, tagname, layer=k, iper=iper) + ! call nf_verify(nf90_def_var(ncid, varname_l, NF90_DOUBLE, & (/dim_ids%nmesh_face/), var_id(k)), & nc_fname) ! ! -- apply chunking parameters - if (chunk_face > 0) then - call nf_verify(nf90_def_var_chunking(ncid, var_id(k), NF90_CHUNKED, & - (/chunk_face/)), nc_fname) - end if + call ncvar_chunk(ncid, var_id(k), chunk_face, nc_fname) ! -- deflate and shuffle - if (deflate >= 0) then - call nf_verify(nf90_def_var_deflate(ncid, var_id(k), shuffle=shuffle, & - deflate=1, deflate_level=deflate), & - nc_fname) - end if + call ncvar_deflate(ncid, var_id(k), deflate, shuffle, nc_fname) ! + ! -- put attr call nf_verify(nf90_put_att(ncid, var_id(k), '_FillValue', & - (/NF90_FILL_DOUBLE/)), nc_fname) + (/fill_value/)), nc_fname) call nf_verify(nf90_put_att(ncid, var_id(k), 'long_name', & longname_l), nc_fname) - if (gridmap_name /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'coordinates', & - 'mesh_face_x mesh_face_y'), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'grid_mapping', & - gridmap_name), nc_fname) - end if - ! - if (nc_tag /= '') then - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_input', & - nc_tag), nc_fname) - call nf_verify(nf90_put_att(ncid, var_id(k), 'modflow6_layer', & - k), nc_fname) - end if + ! + ! -- add grid mapping and mf6 attr + call ncvar_gridmap(ncid, var_id(k), gridmap_name, nc_fname) + call ncvar_mf6attr(ncid, var_id(k), k, iper, iaux, nc_tag, nc_fname) end do ! ! -- exit define mode and write data diff --git a/src/Utilities/Export/MeshNCModel.f90 b/src/Utilities/Export/MeshNCModel.f90 index bb9de03b734..4e4ff3fab09 100644 --- a/src/Utilities/Export/MeshNCModel.f90 +++ b/src/Utilities/Export/MeshNCModel.f90 @@ -22,6 +22,11 @@ module MeshModelModule private public :: MeshNCDimIdType, MeshNCVarIdType public :: Mesh2dModelType + public :: ncvar_chunk + public :: ncvar_deflate + public :: ncvar_gridmap + public :: ncvar_mf6attr + public :: export_varname !> @brief type for storing model export dimension ids !< @@ -242,11 +247,6 @@ subroutine add_pkg_data(this) if (export_arrays > 0) then pkgtype = idm_subcomponent_type(this%modeltype, ptype) param_dfns => param_definitions(this%modeltype, pkgtype) - !if (idm_multi_package(this%modeltype, pkgtype)) then - ! call this%export_input_arrays(pname, mempath, param_dfns) - !else - ! call this%export_input_arrays(pkgtype, mempath, param_dfns) - !end if call this%export_input_arrays(ptype, pname, mempath, param_dfns) end if end if @@ -292,13 +292,8 @@ subroutine define_dependent(this) this%nc_fname) end if ! -- deflate and shuffle - if (this%deflate >= 0) then - call nf_verify(nf90_def_var_deflate(this%ncid, & - this%var_ids%dependent(k), & - shuffle=this%shuffle, deflate=1, & - deflate_level=this%deflate), & - this%nc_fname) - end if + call ncvar_deflate(this%ncid, this%var_ids%dependent(k), this%deflate, & + this%shuffle, this%nc_fname) ! ! -- assign variable attributes call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), & @@ -316,15 +311,9 @@ subroutine define_dependent(this) call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), & 'location', 'face'), this%nc_fname) ! - if (this%ogc_wkt /= '') then - ! -- associate with projection - call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), & - 'coordinates', 'mesh_face_x mesh_face_y'), & - this%nc_fname) - call nf_verify(nf90_put_att(this%ncid, this%var_ids%dependent(k), & - 'grid_mapping', this%gridmap_name), & - this%nc_fname) - end if + ! -- add grid mapping + call ncvar_gridmap(this%ncid, this%var_ids%dependent(k), & + this%gridmap_name, this%nc_fname) end do end subroutine define_dependent @@ -504,4 +493,115 @@ subroutine create_mesh(this) 'start_index', 1), this%nc_fname) end subroutine create_mesh + !> @brief define variable chunking + !< + subroutine ncvar_chunk(ncid, varid, chunk_face, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: chunk_face + character(len=*), intent(in) :: nc_fname + ! + ! -- apply chunking parameters + if (chunk_face > 0) then + call nf_verify(nf90_def_var_chunking(ncid, varid, NF90_CHUNKED, & + (/chunk_face/)), nc_fname) + end if + end subroutine ncvar_chunk + + !> @brief define variable compression + !< + subroutine ncvar_deflate(ncid, varid, deflate, shuffle, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: deflate + integer(I4B), intent(in) :: shuffle + character(len=*), intent(in) :: nc_fname + ! -- deflate and shuffle + if (deflate >= 0) then + call nf_verify(nf90_def_var_deflate(ncid, varid, shuffle=shuffle, & + deflate=1, deflate_level=deflate), & + nc_fname) + end if + end subroutine ncvar_deflate + + !> @brief put variable gridmap attributes + !< + subroutine ncvar_gridmap(ncid, varid, gridmap_name, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + character(len=*), intent(in) :: gridmap_name + character(len=*), intent(in) :: nc_fname + ! + if (gridmap_name /= '') then + call nf_verify(nf90_put_att(ncid, varid, 'coordinates', & + 'mesh_face_x mesh_face_y'), nc_fname) + call nf_verify(nf90_put_att(ncid, varid, 'grid_mapping', & + gridmap_name), nc_fname) + end if + end subroutine ncvar_gridmap + + !> @brief put variable internal attributes + !< + subroutine ncvar_mf6attr(ncid, varid, layer, iper, iaux, nc_tag, nc_fname) + integer(I4B), intent(in) :: ncid + integer(I4B), intent(in) :: varid + integer(I4B), intent(in) :: layer + integer(I4B), intent(in) :: iper + integer(I4B), intent(in) :: iaux + character(len=*), intent(in) :: nc_tag + character(len=*), intent(in) :: nc_fname + ! + if (nc_tag /= '') then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_input', & + nc_tag), nc_fname) + if (layer > 0) then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_layer', & + layer), nc_fname) + end if + ! + if (iper > 0) then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_iper', & + iper), nc_fname) + end if + ! + if (iaux > 0) then + call nf_verify(nf90_put_att(ncid, varid, 'modflow6_iaux', & + iaux), nc_fname) + end if + end if + end subroutine ncvar_mf6attr + + !> @brief build netcdf variable name + !< + function export_varname(varname, layer, iper, iaux) result(vname) + use InputOutputModule, only: lowcase + character(len=*), intent(in) :: varname + integer(I4B), optional, intent(in) :: layer + integer(I4B), optional, intent(in) :: iper + integer(I4B), optional, intent(in) :: iaux + character(len=LINELENGTH) :: vname + ! + vname = '' + ! + if (varname /= '') then + vname = varname + call lowcase(vname) + if (present(layer)) then + if (layer > 0) then + write (vname, '(a,i0)') trim(vname)//'_l', layer + end if + end if + if (present(iper)) then + if (iper > 0) then + write (vname, '(a,i0)') trim(vname)//'_p', iper + end if + end if + if (present(iaux)) then + if (iaux > 0) then + write (vname, '(a,i0)') trim(vname)//'a', iaux + end if + end if + end if + end function export_varname + end module MeshModelModule diff --git a/src/Utilities/Export/ModelExport.f90 b/src/Utilities/Export/ModelExport.f90 index cf01bc28f51..74a60fb4771 100644 --- a/src/Utilities/Export/ModelExport.f90 +++ b/src/Utilities/Export/ModelExport.f90 @@ -12,10 +12,12 @@ module ModelExportModule LENMEMPATH use ListModule, only: ListType use NCModelExportModule, only: NCBaseModelExportType + use InputLoadTypeModule, only: ModelDynamicPkgsType implicit none private public :: modelexports_create + public :: modelexports_post_prepare public :: modelexports_post_step public :: modelexports_destroy public :: nc_export_active @@ -32,6 +34,7 @@ module ModelExportModule !! !< type :: ExportModelType + type(ModelDynamicPkgsType), pointer :: loaders => null() character(len=LENMODELNAME) :: modelname !< name of model character(len=LENCOMPONENTNAME) :: modeltype !< type of model character(len=LINELENGTH) :: modelfname !< name of model input file @@ -41,6 +44,7 @@ module ModelExportModule integer(I4B) :: iout !< lst file descriptor contains procedure :: init + procedure :: post_prepare procedure :: post_step procedure :: destroy end type ExportModelType @@ -68,7 +72,6 @@ end function nc_export_active !> @brief create export container variable for all local models !! subroutine modelexports_create(iout) - use InputLoadTypeModule, only: ModelDynamicPkgsType use InputLoadTypeModule, only: model_dynamic_pkgs use MemoryManagerModule, only: mem_setptr use MemoryManagerExtModule, only: mem_set_value @@ -103,9 +106,7 @@ subroutine modelexports_create(iout) call mem_setptr(disenum, 'DISENUM', model_mempath) ! ! -- initialize model - call export_model%init(model_dynamic_input%modelname, & - model_dynamic_input%modeltype, & - model_dynamic_input%modelfname, disenum, iout) + call export_model%init(model_dynamic_input, disenum, iout) ! ! -- update EXPORT_NETCDF string if provided call mem_set_value(exportstr, 'EXPORT_NETCDF', modelnam_mempath, found) @@ -113,6 +114,7 @@ subroutine modelexports_create(iout) if (exportstr == 'STRUCTURED') then export_model%nctype = NETCDF_STRUCTURED else + ! -- mesh export is default export_model%nctype = NETCDF_UGRID end if end if @@ -122,6 +124,26 @@ subroutine modelexports_create(iout) end do end subroutine modelexports_create + !> @brief export model list post prepare step + !! + subroutine modelexports_post_prepare() + ! -- local variables + class(*), pointer :: obj + class(ExportModelType), pointer :: export_model + integer(I4B) :: n + ! + do n = 1, export_models%Count() + obj => export_models%GetItem(n) + if (associated(obj)) then + select type (obj) + class is (ExportModelType) + export_model => obj + call export_model%post_prepare() + end select + end if + end do + end subroutine modelexports_post_prepare + !> @brief export model list post step !! subroutine modelexports_post_step() @@ -169,18 +191,17 @@ end subroutine modelexports_destroy !> @brief initialize model export container variable !! !< - subroutine init(this, modelname, modeltype, modelfname, disenum, iout) + subroutine init(this, loaders, disenum, iout) use NCModelExportModule, only: NETCDF_UNDEF class(ExportModelType), intent(inout) :: this - character(len=*), intent(in) :: modelname - character(len=*), intent(in) :: modeltype - character(len=*), intent(in) :: modelfname + type(ModelDynamicPkgsType), pointer, intent(in) :: loaders integer(I4B), intent(in) :: disenum integer(I4B), intent(in) :: iout ! - this%modelname = modelname - this%modeltype = modeltype - this%modelfname = modelfname + this%loaders => loaders + this%modelname = loaders%modelname + this%modeltype = loaders%modeltype + this%modelfname = loaders%modelfname this%nctype = NETCDF_UNDEF this%disenum = disenum this%iout = iout @@ -188,6 +209,17 @@ subroutine init(this, modelname, modeltype, modelfname, disenum, iout) nullify (this%nc_export) end subroutine init + !> @brief model export container post prepare step actions + !! + !< + subroutine post_prepare(this) + class(ExportModelType), intent(inout) :: this + ! + if (associated(this%nc_export)) then + call this%nc_export%export_input() + end if + end subroutine post_prepare + !> @brief model export container post step actions !! !< diff --git a/src/Utilities/Export/NCExportCreate.f90 b/src/Utilities/Export/NCExportCreate.f90 index 3ce0a33970d..c10c74bb3a2 100644 --- a/src/Utilities/Export/NCExportCreate.f90 +++ b/src/Utilities/Export/NCExportCreate.f90 @@ -17,6 +17,7 @@ module NCExportCreateModule use DisuModule, only: DisuType use ModelExportModule, only: export_models, get_export_model use ModelExportModule, only: ExportModelType + use NCModelExportModule, only: ExportPackageType implicit none private @@ -31,6 +32,7 @@ subroutine create_nc_export(export_model, num_model) use MeshDisModelModule, only: Mesh2dDisExportType use MeshDisvModelModule, only: Mesh2dDisvExportType use DisNCStructuredModule, only: DisNCStructuredType + use InputLoadTypeModule, only: ModelDynamicPkgsType type(ExportModelType), pointer, intent(inout) :: export_model class(NumericalModelType), pointer, intent(in) :: num_model class(Mesh2dDisExportType), pointer :: ugrid_dis @@ -53,6 +55,10 @@ subroutine create_nc_export(export_model, num_model) ugrid_dis%dis => disbase end select ! + ! -- set dynamic loaders + call create_export_pkglist(ugrid_dis%pkglist, export_model%loaders, & + export_model%iout) + ! ! -- initialize export object call ugrid_dis%init(export_model%modelname, export_model%modeltype, & export_model%modelfname, export_model%disenum, & @@ -75,6 +81,10 @@ subroutine create_nc_export(export_model, num_model) structured_dis%dis => disbase end select ! + ! -- set dynamic loaders + call create_export_pkglist(structured_dis%pkglist, export_model%loaders, & + export_model%iout) + ! ! -- initialize export object call structured_dis%init(export_model%modelname, export_model%modeltype, & export_model%modelfname, export_model%disenum, & @@ -98,6 +108,10 @@ subroutine create_nc_export(export_model, num_model) ugrid_disv%disv => disbase end select ! + ! -- set dynamic loaders + call create_export_pkglist(ugrid_disv%pkglist, export_model%loaders, & + export_model%iout) + ! ! -- initialize export object call ugrid_disv%init(export_model%modelname, export_model%modeltype, & export_model%modelfname, export_model%disenum, & @@ -123,6 +137,64 @@ subroutine create_nc_export(export_model, num_model) end select end subroutine create_nc_export + subroutine create_export_pkglist(pkglist, loaders, iout) + use ListModule, only: ListType + use MemoryManagerExtModule, only: mem_set_value + use InputLoadTypeModule, only: ModelDynamicPkgsType + use InputLoadTypeModule, only: DynamicPkgLoadBaseType + use AsciiInputLoadTypeModule, only: AsciiDynamicPkgLoadBaseType + use Mf6FileGridInputModule, only: BoundGridInputType + use IdmMf6FileModule, only: Mf6FileDynamicPkgLoadType + type(ListType), intent(inout) :: pkglist + type(ModelDynamicPkgsType), pointer, intent(in) :: loaders + integer(I4B), intent(in) :: iout + class(DynamicPkgLoadBaseType), pointer :: dynamic_pkg + class(AsciiDynamicPkgLoadBaseType), pointer :: rp_loader + type(ExportPackageType), pointer :: export_pkg + integer(I4B), pointer :: export_arrays + class(*), pointer :: obj + logical(LGP) :: found + integer(I4B) :: n + ! + ! -- create list of in scope loaders + allocate (export_arrays) + ! + do n = 1, loaders%pkglist%Count() + ! + ! -- initialize export arrays option + export_arrays = 0 + ! + dynamic_pkg => loaders%get(n) + ! + ! -- update export arrays option + call mem_set_value(export_arrays, 'EXPORT_NC', & + dynamic_pkg%mf6_input%mempath, found) + ! + if (export_arrays > 0 .and. dynamic_pkg%readasarrays) then + select type (dynamic_pkg) + type is (Mf6FileDynamicPkgLoadType) + ! + rp_loader => dynamic_pkg%rp_loader + ! + select type (rp_loader) + type is (BoundGridInputType) + ! -- create the export object + allocate (export_pkg) + call export_pkg%init(rp_loader%mf6_input, & + rp_loader%bound_context%mshape, & + rp_loader%param_names, rp_loader%nparam) + obj => export_pkg + call pkglist%add(obj) + ! + end select + end select + end if + end do + ! + ! -- cleanup + deallocate (export_arrays) + end subroutine create_export_pkglist + !> @brief initialize netcdf model export type !! subroutine nc_export_create() diff --git a/src/Utilities/Export/NCModel.f90 b/src/Utilities/Export/NCModel.f90 index 58fb8d467e0..cb6435e557f 100644 --- a/src/Utilities/Export/NCModel.f90 +++ b/src/Utilities/Export/NCModel.f90 @@ -9,24 +9,43 @@ module NCModelExportModule use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENCOMPONENTNAME, LENMODELNAME, & - LENMEMPATH, LENBIGLINE, DIS, DISU, DISV + LENMEMPATH, LENBIGLINE, LENVARNAME, DIS, DISU, DISV use SimVariablesModule, only: errmsg use SimModule, only: store_error, store_error_filename + use InputLoadTypeModule, only: ModelDynamicPkgsType + use ModflowInputModule, only: ModflowInputType + use BoundInputContextModule, only: ReadStateVarType, rsv_name + use ListModule, only: ListType implicit none private public :: NCBaseModelExportType, NCModelExportType public :: NCExportAnnotation + public :: ExportPackageType public :: NETCDF_UNDEF, NETCDF_STRUCTURED, NETCDF_UGRID + public :: export_longname !> @brief netcdf export types enumerator !< ENUM, BIND(C) ENUMERATOR :: NETCDF_UNDEF = 0 !< undefined netcdf export type - ENUMERATOR :: NETCDF_STRUCTURED = 1 !< netcdf structrured export - ENUMERATOR :: NETCDF_UGRID = 2 !< netcdf mesh export + ENUMERATOR :: NETCDF_UGRID = 1 !< netcdf mesh export + ENUMERATOR :: NETCDF_STRUCTURED = 2 !< netcdf structrured export END ENUM + type :: ExportPackageType + type(ModflowInputType) :: mf6_input !< description of modflow6 input + character(len=LINELENGTH), dimension(:), allocatable :: param_names !< dynamic param tagnames + type(ReadStateVarType), dimension(:), allocatable :: param_reads !< param read states + integer(I4B), dimension(:), pointer, contiguous :: mshape => null() !< model shape + integer(I4B), pointer :: iper !< most recent package rp load + integer(I4B) :: iper_export !< most recent period of netcdf package export + integer(I4B) :: nparam !< number of in scope params + contains + procedure :: init => epkg_init + procedure :: destroy => epkg_destroy + end type ExportPackageType + !> @brief netcdf export attribute annotations !< type :: NCExportAnnotation @@ -45,11 +64,12 @@ module NCModelExportModule !> @brief base class for an export model !< type :: NCModelExportType + type(ListType) :: pkglist character(len=LENMODELNAME) :: modelname !< name of model character(len=LENCOMPONENTNAME) :: modeltype !< type of model character(len=LINELENGTH) :: modelfname !< name of model input file character(len=LINELENGTH) :: nc_fname !< name of netcdf export file - character(len=LINELENGTH) :: gridmap_name = 'projection' !< name of grid mapping variable + character(len=LINELENGTH) :: gridmap_name !< name of grid mapping variable character(len=LINELENGTH) :: mesh_name = 'mesh' !< name of mesh container variable character(len=LENMEMPATH) :: dis_mempath !< discretization input mempath character(len=LENMEMPATH) :: ncf_mempath !< netcdf utility package input mempath @@ -70,6 +90,8 @@ module NCModelExportModule logical(LGP) :: chunking_active !< have chunking parameters been provided contains procedure :: init => export_init + procedure :: get => export_get + procedure :: input_attribute procedure :: destroy => export_destroy end type NCModelExportType @@ -77,25 +99,95 @@ module NCModelExportModule !< type, abstract, extends(NCModelExportType) :: NCBaseModelExportType contains - procedure(model_define_if), deferred :: df - procedure(model_step_if), deferred :: step + procedure :: export_input + procedure(model_define), deferred :: df + procedure(model_step), deferred :: step + procedure(package_export), deferred :: package_step + procedure(package_export_ilayer), deferred :: package_step_ilayer end type NCBaseModelExportType !> @brief abstract interfaces for model netcdf export type !< abstract interface - subroutine model_define_if(this) + subroutine model_define(this) import NCBaseModelExportType class(NCBaseModelExportType), intent(inout) :: this end subroutine - subroutine model_step_if(this) + subroutine model_step(this) import NCBaseModelExportType class(NCBaseModelExportType), intent(inout) :: this end subroutine + subroutine package_export(this, export_pkg) + import NCBaseModelExportType, ExportPackageType + class(NCBaseModelExportType), intent(inout) :: this + class(ExportPackageType), intent(in) :: export_pkg + end subroutine + subroutine package_export_ilayer(this, export_pkg, ilayer_varname, & + ilayer) + import NCBaseModelExportType, ExportPackageType, I4B + class(NCBaseModelExportType), intent(inout) :: this + class(ExportPackageType), intent(in) :: export_pkg + character(len=*), intent(in) :: ilayer_varname + integer(I4B), intent(in) :: ilayer + end subroutine end interface contains + !> @brief initialize dynamic package export object + !< + subroutine epkg_init(this, mf6_input, mshape, param_names, & + nparam) + use SimVariablesModule, only: idm_context + use MemoryManagerModule, only: mem_setptr + use MemoryManagerExtModule, only: mem_set_value + use MemoryHelperModule, only: create_mem_path + ! -- dummy + class(ExportPackageType), intent(inout) :: this + type(ModflowInputType), intent(in) :: mf6_input + integer(I4B), dimension(:), pointer, contiguous, intent(in) :: mshape !< model shape + character(len=LINELENGTH), dimension(:), allocatable, & + intent(in) :: param_names + integer(I4B), intent(in) :: nparam + integer(I4B) :: n + character(len=LENVARNAME) :: rs_varname + character(len=LENMEMPATH) :: input_mempath + integer(I4B), pointer :: rsvar + ! + this%mf6_input = mf6_input + this%mshape => mshape + this%nparam = nparam + this%iper_export = 0 + ! + input_mempath = create_mem_path(component=mf6_input%component_name, & + subcomponent=mf6_input%subcomponent_name, & + context=idm_context) + ! + ! -- allocate param arrays + allocate (this%param_names(nparam)) + allocate (this%param_reads(nparam)) + ! + ! -- set param arrays + do n = 1, nparam + this%param_names(n) = param_names(n) + rs_varname = rsv_name(param_names(n)) + call mem_setptr(rsvar, rs_varname, mf6_input%mempath) + this%param_reads(n)%invar => rsvar + end do + ! + ! -- set pointer to loaded input period + call mem_setptr(this%iper, 'IPER', mf6_input%mempath) + end subroutine epkg_init + + !> @brief destroy dynamic package export object + !< + subroutine epkg_destroy(this) + use InputDefinitionModule, only: InputParamDefinitionType + ! -- dummy + class(ExportPackageType), intent(inout) :: this + if (allocated(this%param_names)) deallocate (this%param_names) + end subroutine epkg_destroy + !> @brief set netcdf file scoped attributes !< subroutine set(this, modelname, modeltype, modelfname, nctype) @@ -197,7 +289,7 @@ subroutine export_init(this, modelname, modeltype, modelfname, disenum, & this%modeltype = modeltype this%modelfname = modelfname this%nc_fname = trim(modelname)//'.nc' - call lowcase(this%nc_fname) + this%gridmap_name = '' this%ncf_mempath = '' this%ogc_wkt = '' this%datetime = '' @@ -213,6 +305,8 @@ subroutine export_init(this, modelname, modeltype, modelfname, disenum, & this%iout = iout this%chunking_active = .false. ! + call lowcase(this%nc_fname) + ! ! -- set file scoped attributes call this%annotation%set(modelname, modeltype, modelfname, nctype) ! @@ -260,6 +354,10 @@ subroutine export_init(this, modelname, modeltype, modelfname, disenum, & found%chunk_time) end if ! + if (found%ogc_wkt) then + this%gridmap_name = 'projection' + end if + ! ! -- ATTR_OFF turns off modflow 6 input attributes if (found%attr_off) then this%input_attr = 0 @@ -278,12 +376,121 @@ subroutine export_init(this, modelname, modeltype, modelfname, disenum, & this%totnstp = sum(nstp) end subroutine export_init + !> @brief retrieve dynamic export object from package list + !< + function export_get(this, idx) result(res) + use ListModule, only: ListType + class(NCModelExportType), intent(inout) :: this + integer(I4B), intent(in) :: idx + class(ExportPackageType), pointer :: res + class(*), pointer :: obj + ! + nullify (res) + obj => this%pkglist%GetItem(idx) + if (associated(obj)) then + select type (obj) + class is (ExportPackageType) + res => obj + end select + end if + end function export_get + + !> @brief build modflow6_input attribute string + !< + function input_attribute(this, pkgname, idt) result(attr) + use InputOutputModule, only: lowcase + use MemoryHelperModule, only: memPathSeparator + use InputDefinitionModule, only: InputParamDefinitionType + class(NCModelExportType), intent(inout) :: this + character(len=*), intent(in) :: pkgname + type(InputParamDefinitionType), pointer, intent(in) :: idt + character(len=LINELENGTH) :: attr + ! + attr = '' + ! + if (this%input_attr > 0) then + attr = trim(this%modelname)//memPathSeparator//trim(pkgname)// & + memPathSeparator//trim(idt%mf6varname) + end if + end function input_attribute + + !> @brief build netcdf variable longname + !< + function export_longname(longname, pkgname, tagname, layer, iper) result(lname) + use InputOutputModule, only: lowcase + character(len=*), intent(in) :: longname + character(len=*), intent(in) :: pkgname + character(len=*), intent(in) :: tagname + integer(I4B), intent(in) :: layer + integer(I4B), optional, intent(in) :: iper + character(len=LINELENGTH) :: lname + character(len=LINELENGTH) :: pname, vname + ! + pname = pkgname + vname = tagname + call lowcase(pname) + call lowcase(vname) + if (longname == '') then + lname = trim(pname)//' '//trim(vname) + else + lname = longname + end if + if (layer > 0) then + write (lname, '(a,i0)') trim(lname)//' layer=', layer + end if + if (present(iper)) then + if (iper > 0) then + write (lname, '(a,i0)') trim(lname)//' period=', iper + end if + end if + end function export_longname + + !> @brief netcdf dynamic package period export + !< + subroutine export_input(this) + use TdisModule, only: kper + use ArrayHandlersModule, only: ifind + class(NCBaseModelExportType), intent(inout) :: this + integer(I4B) :: idx, ilayer + class(ExportPackageType), pointer :: export_pkg + character(len=LENVARNAME) :: ilayer_varname + ! + do idx = 1, this%pkglist%Count() + ! + export_pkg => this%get(idx) + ! -- last loaded data is not current period + if (export_pkg%iper /= kper) cycle + ! -- period input already exported + if (export_pkg%iper_export >= export_pkg%iper) cycle + ! -- set exported iper + export_pkg%iper_export = export_pkg%iper + ! + ! -- initialize ilayer + ilayer = 0 + ! + ! -- set expected ilayer index variable name + ilayer_varname = 'I'//trim(export_pkg%mf6_input%subcomponent_type(1:3)) + ! + ! -- is ilayer variable in param name list + ilayer = ifind(export_pkg%param_names, ilayer_varname) + ! + ! -- layer index variable is required to be first defined in period block + if (ilayer == 1) then + call this%package_step_ilayer(export_pkg, ilayer_varname, ilayer) + else + call this%package_step(export_pkg) + end if + ! + end do + end subroutine export_input + !> @brief destroy model netcdf export object !< subroutine export_destroy(this) use MemoryManagerExtModule, only: memorystore_remove use SimVariablesModule, only: idm_context class(NCModelExportType), intent(inout) :: this + ! ! -- override in derived class deallocate (this%deflate) deallocate (this%shuffle) diff --git a/src/Utilities/Idm/BoundInputContext.f90 b/src/Utilities/Idm/BoundInputContext.f90 index cecc3d7282d..bf2e42c5948 100644 --- a/src/Utilities/Idm/BoundInputContext.f90 +++ b/src/Utilities/Idm/BoundInputContext.f90 @@ -20,6 +20,7 @@ module BoundInputContextModule private public :: BoundInputContextType public :: ReadStateVarType + public :: rsv_name !> @brief Pointer type for read state variable !< @@ -357,18 +358,9 @@ function rsv_alloc(this, mf6varname) result(varname) character(len=*), intent(in) :: mf6varname ! -- local character(len=LENVARNAME) :: varname - integer(I4B) :: ilen integer(I4B), pointer :: intvar - character(len=2) :: prefix = 'IN' - ! - ! -- assign first column as the block number - ilen = len_trim(mf6varname) ! - if (ilen > (LENVARNAME - len(prefix))) then - varname = prefix//mf6varname(1:(LENVARNAME - len(prefix))) - else - varname = prefix//trim(mf6varname) - end if + varname = rsv_name(mf6varname) ! call mem_allocate(intvar, varname, this%mf6_input%mempath) intvar = -1 @@ -424,4 +416,29 @@ subroutine bound_params(this, params, nparam, input_name, create) return end subroutine bound_params + !> @brief create read state variable name + !! + !< + function rsv_name(mf6varname) result(varname) + ! -- modules + use ConstantsModule, only: LENVARNAME + ! -- dummy + character(len=*), intent(in) :: mf6varname + ! -- local + character(len=LENVARNAME) :: varname + integer(I4B) :: ilen + character(len=2) :: prefix = 'IN' + ! + ilen = len_trim(mf6varname) + ! + if (ilen > (LENVARNAME - len(prefix))) then + varname = prefix//mf6varname(1:(LENVARNAME - len(prefix))) + else + varname = prefix//trim(mf6varname) + end if + ! + ! -- return + return + end function rsv_name + end module BoundInputContextModule diff --git a/src/Utilities/Idm/IdmLoad.f90 b/src/Utilities/Idm/IdmLoad.f90 index 654c60f0232..4824a47d3d8 100644 --- a/src/Utilities/Idm/IdmLoad.f90 +++ b/src/Utilities/Idm/IdmLoad.f90 @@ -583,7 +583,7 @@ subroutine dynamic_da(iout) return end subroutine dynamic_da - !> @brief return sim input context PRINT_INTPUT value + !> @brief return sim input context PRINT_INPUT value !< function input_param_log() result(paramlog) use MemoryHelperModule, only: create_mem_path diff --git a/src/Utilities/Idm/ModelPackageInputs.f90 b/src/Utilities/Idm/ModelPackageInputs.f90 index 8ca93fa584d..24f35cdb88e 100644 --- a/src/Utilities/Idm/ModelPackageInputs.f90 +++ b/src/Utilities/Idm/ModelPackageInputs.f90 @@ -181,7 +181,13 @@ subroutine pkgtype_add(this, modelname, mtype_component, filetype, & ! ! -- set pkgname if empty if (this%pkgnames(this%pnum) == '') then - write (pname, '(a,i0)') trim(this%subcomponent_type)//'-', this%pnum + if (multi_pkg_type(mtype_component, & + this%subcomponent_type, & + filetype)) then + write (pname, '(a,i0)') trim(this%subcomponent_type)//'-', this%pnum + else + write (pname, '(a,i0)') trim(this%subcomponent_type) + end if this%pkgnames(this%pnum) = pname end if ! diff --git a/src/Utilities/Idm/SourceLoad.F90 b/src/Utilities/Idm/SourceLoad.F90 index e8f7f55c3ce..3090e42232d 100644 --- a/src/Utilities/Idm/SourceLoad.F90 +++ b/src/Utilities/Idm/SourceLoad.F90 @@ -21,7 +21,8 @@ module SourceLoadModule public :: open_source_file public :: load_modelnam, load_simnam, load_simtdis public :: remote_model_ndim - public :: export_cr, export_post_step, export_da + public :: export_cr, export_da + public :: export_post_prepare, export_post_step public :: nc_close public :: netcdf_context @@ -65,7 +66,7 @@ function create_input_loader(component_type, subcomponent_type, & ! -- initialize loader call loader%init(mf6_input, component_name, component_fname, input_fname) ! - ! -- initialize loader netcdf variables structure + ! -- initialize loader netcdf variables data structure if (present(nc_vars)) then call nc_vars%create_varlists(component_name, sc_name, loader%nc_vars) else @@ -362,6 +363,14 @@ subroutine export_cr() end if end subroutine export_cr + !> @brief model exports post prepare step actions + !< + subroutine export_post_prepare() + ! -- modules + use ModelExportModule, only: modelexports_post_prepare + call modelexports_post_prepare() + end subroutine export_post_prepare + !> @brief model exports post step actions !< subroutine export_post_step() @@ -406,7 +415,7 @@ function netcdf_context(modeltype, component_type, modelname, & #if defined(__WITH_NETCDF__) use NCContextBuildModule, only: open_ncfile, create_netcdf_context #endif - ! -- drummy + ! -- dummy character(len=*), intent(in) :: modeltype character(len=*), intent(in) :: component_type character(len=*), intent(in) :: modelname diff --git a/src/Utilities/Idm/mf6blockfile/LoadNCInput.F90 b/src/Utilities/Idm/mf6blockfile/LoadNCInput.F90 index ef9c9986c81..f2164d3e421 100644 --- a/src/Utilities/Idm/mf6blockfile/LoadNCInput.F90 +++ b/src/Utilities/Idm/mf6blockfile/LoadNCInput.F90 @@ -45,7 +45,7 @@ end subroutine error_and_exit !> @brief Read a NetCDF integer 1d array !< subroutine nc_read_int1d(int1d, mshape, idt, mf6_input, nc_vars, input_fname, & - iout) + iout, kper) integer(I4B), dimension(:), pointer, contiguous, intent(in) :: int1d integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape type(InputParamDefinitionType), intent(in) :: idt !< input data type object describing this record @@ -53,9 +53,10 @@ subroutine nc_read_int1d(int1d, mshape, idt, mf6_input, nc_vars, input_fname, & type(NCPackageVarsType), pointer, intent(in) :: nc_vars character(len=*), intent(in) :: input_fname integer(I4B), intent(in) :: iout + integer(I4B), optional, intent(in) :: kper #if defined(__WITH_NETCDF__) call netcdf_array_load(int1d, mshape, idt, mf6_input, nc_vars, input_fname, & - iout) + iout, kper) #else call error_and_exit(idt%tagname, input_fname) #endif @@ -102,7 +103,7 @@ end subroutine nc_read_int3d !> @brief Read a NetCDF double 1d array !< subroutine nc_read_dbl1d(dbl1d, mshape, idt, mf6_input, nc_vars, input_fname, & - iout) + iout, kper, iaux) real(DP), dimension(:), pointer, contiguous, intent(in) :: dbl1d integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape type(InputParamDefinitionType), intent(in) :: idt !< input data type object describing this record @@ -110,9 +111,11 @@ subroutine nc_read_dbl1d(dbl1d, mshape, idt, mf6_input, nc_vars, input_fname, & type(NCPackageVarsType), pointer, intent(in) :: nc_vars character(len=*), intent(in) :: input_fname integer(I4B), intent(in) :: iout + integer(I4B), optional, intent(in) :: kper + integer(I4B), optional, intent(in) :: iaux #if defined(__WITH_NETCDF__) call netcdf_array_load(dbl1d, mshape, idt, mf6_input, nc_vars, input_fname, & - iout) + iout, kper=kper, iaux=iaux) #else call error_and_exit(idt%tagname, input_fname) #endif diff --git a/src/Utilities/Idm/mf6blockfile/Mf6FileGridInput.f90 b/src/Utilities/Idm/mf6blockfile/Mf6FileGridInput.f90 index 3fcb375b1ef..6e87153da4f 100644 --- a/src/Utilities/Idm/mf6blockfile/Mf6FileGridInput.f90 +++ b/src/Utilities/Idm/mf6blockfile/Mf6FileGridInput.f90 @@ -153,11 +153,12 @@ subroutine bndgrid_rp(this, parser) class(BoundGridInputType), intent(inout) :: this !< Mf6FileGridInputType type(BlockParserType), pointer, intent(inout) :: parser ! -- local - logical(LGP) :: endOfBlock + logical(LGP) :: endOfBlock, netcdf character(len=LINELENGTH) :: keyword, param_tag type(InputParamDefinitionType), pointer :: idt integer(I4B) :: iaux, iparam character(len=LENTIMESERIESNAME) :: tas_name + integer(I4B), dimension(:), pointer, contiguous :: int1d ! ! -- reset for this period call this%reset() @@ -170,6 +171,7 @@ subroutine bndgrid_rp(this, parser) do ! -- initialize iaux = 0 + netcdf = .false. ! ! -- read next line call parser%GetNextLine(endOfBlock) @@ -190,11 +192,12 @@ subroutine bndgrid_rp(this, parser) this%mf6_input%subcomponent_type, & 'PERIOD', param_tag, this%input_name) ! - ! -- look for TAS keyword if tas is active - if (this%tas_active /= 0) then - call parser%GetStringCaps(keyword) - ! - if (keyword == 'TIMEARRAYSERIES') then + ! -- look for TAS and NetCDF keywords + call parser%GetStringCaps(keyword) + ! + if (keyword == 'TIMEARRAYSERIES') then + if (this%tas_active /= 0) then + ! call parser%GetStringCaps(tas_name) ! if (param_tag == 'AUX') then @@ -210,16 +213,32 @@ subroutine bndgrid_rp(this, parser) ! ! -- cycle to next input param cycle + else + ! TODO: throw error end if ! + else if (keyword == 'NETCDF') then + netcdf = .true. end if ! ! -- read and load the parameter - call this%param_load(parser, idt%datatype, idt%mf6varname, idt%tagname, & - idt%shape, this%mf6_input%mempath, iaux) + call this%param_load(parser, idt, this%mf6_input%mempath, netcdf, iaux) ! end do ! + ! -- check if layer index variable was read + ! TODO: assumes layer index variable is always in scope + if (this%param_reads(1)%invar == 0) then + ! -- set to default of 1 without updating invar + idt => get_param_definition_type(this%mf6_input%param_dfns, & + this%mf6_input%component_type, & + this%mf6_input%subcomponent_type, & + 'PERIOD', this%param_names(1), & + this%input_name) + call mem_setptr(int1d, idt%mf6varname, this%mf6_input%mempath) + int1d = 1 + end if + ! ! if (this%tas_active /= 0) then call this%tas_links_create(parser%iuactive) @@ -322,9 +341,9 @@ subroutine bndgrid_params_alloc(this) return end subroutine bndgrid_params_alloc - subroutine bndgrid_param_load(this, parser, datatype, varname, & - tagname, shapestr, mempath, iaux) + subroutine bndgrid_param_load(this, parser, idt, mempath, netcdf, iaux) ! -- modules + use TdisModule, only: kper use MemoryManagerModule, only: mem_setptr use ArrayHandlersModule, only: ifind use InputDefinitionModule, only: InputParamDefinitionType @@ -332,15 +351,14 @@ subroutine bndgrid_param_load(this, parser, datatype, varname, & use Double1dReaderModule, only: read_dbl1d use Double2dReaderModule, only: read_dbl2d use Integer1dReaderModule, only: read_int1d + use LoadNCInputModule, only: netcdf_read_array use IdmLoggerModule, only: idm_log_var ! -- dummy class(BoundGridInputType), intent(inout) :: this !< BoundGridInputType type(BlockParserType), intent(in) :: parser - character(len=*), intent(in) :: datatype - character(len=*), intent(in) :: varname - character(len=*), intent(in) :: tagname - character(len=*), intent(in) :: shapestr + type(InputParamDefinitionType), intent(in) :: idt character(len=*), intent(in) :: mempath + logical(LGP), intent(in) :: netcdf integer(I4B), intent(in) :: iaux ! -- local integer(I4B), dimension(:), pointer, contiguous :: int1d @@ -348,42 +366,59 @@ subroutine bndgrid_param_load(this, parser, datatype, varname, & real(DP), dimension(:, :), pointer, contiguous :: dbl2d integer(I4B) :: iparam, n ! - select case (datatype) + select case (idt%datatype) case ('INTEGER1D') ! - call mem_setptr(int1d, varname, mempath) - call read_int1d(parser, int1d, varname) - call idm_log_var(int1d, tagname, mempath, this%iout) + call mem_setptr(int1d, idt%mf6varname, mempath) + if (netcdf) then + call netcdf_read_array(int1d, this%bound_context%mshape, idt, & + this%mf6_input, this%nc_vars, this%input_name, & + this%iout, kper) + else + call read_int1d(parser, int1d, idt%mf6varname) + end if + call idm_log_var(int1d, idt%tagname, mempath, this%iout) ! case ('DOUBLE1D') ! - call mem_setptr(dbl1d, varname, mempath) - call read_dbl1d(parser, dbl1d, varname) - call idm_log_var(dbl1d, tagname, mempath, this%iout) + call mem_setptr(dbl1d, idt%mf6varname, mempath) + if (netcdf) then + call netcdf_read_array(dbl1d, this%bound_context%mshape, idt, & + this%mf6_input, this%nc_vars, this%input_name, & + this%iout, kper) + else + call read_dbl1d(parser, dbl1d, idt%mf6varname) + end if + call idm_log_var(dbl1d, idt%tagname, mempath, this%iout) ! case ('DOUBLE2D') ! - call mem_setptr(dbl2d, varname, mempath) + call mem_setptr(dbl2d, idt%mf6varname, mempath) allocate (dbl1d(this%bound_context%ncpl)) - call read_dbl1d(parser, dbl1d, varname) + if (netcdf) then + call netcdf_read_array(dbl1d, this%bound_context%mshape, idt, & + this%mf6_input, this%nc_vars, this%input_name, & + this%iout, kper, iaux) + else + call read_dbl1d(parser, dbl1d, idt%mf6varname) + end if do n = 1, this%bound_context%ncpl dbl2d(iaux, n) = dbl1d(n) end do - call idm_log_var(dbl1d, tagname, mempath, this%iout) + call idm_log_var(dbl1d, idt%tagname, mempath, this%iout) deallocate (dbl1d) ! case default ! errmsg = 'IDM unimplemented. Mf6FileGridInput::param_load & - &datatype='//trim(datatype) + &datatype='//trim(idt%datatype) call store_error(errmsg) call store_error_filename(this%input_name) ! end select ! - iparam = ifind(this%param_names, varname) - ! ! -- if param is tracked set read state + iparam = ifind(this%param_names, idt%tagname) if (iparam > 0) then this%param_reads(iparam)%invar = 1 end if diff --git a/src/Utilities/Idm/netcdf/NCArrayReader.f90 b/src/Utilities/Idm/netcdf/NCArrayReader.f90 index e12e76c6843..f728e479a5e 100644 --- a/src/Utilities/Idm/netcdf/NCArrayReader.f90 +++ b/src/Utilities/Idm/netcdf/NCArrayReader.f90 @@ -52,7 +52,7 @@ end function is_layered !> @brief Load NetCDF integer 1D array !< subroutine nc_array_load_int1d(int1d, mshape, idt, mf6_input, nc_vars, & - input_fname, iout) + input_fname, iout, kper) integer(I4B), dimension(:), pointer, contiguous, intent(in) :: int1d integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape type(InputParamDefinitionType), intent(in) :: idt !< input data type object describing this record @@ -60,6 +60,7 @@ subroutine nc_array_load_int1d(int1d, mshape, idt, mf6_input, nc_vars, & type(NCPackageVarsType), pointer, intent(in) :: nc_vars character(len=*), intent(in) :: input_fname integer(I4B), intent(in) :: iout + integer(I4B), optional, intent(in) :: kper ! -- local integer(I4B) :: varid logical(LGP) :: layered @@ -70,7 +71,11 @@ subroutine nc_array_load_int1d(int1d, mshape, idt, mf6_input, nc_vars, & call load_integer1d_layered(int1d, mf6_input, mshape, idt, nc_vars, & input_fname) else - varid = nc_vars%varid(idt%tagname) + if (present(kper)) then + varid = nc_vars%varid(idt%mf6varname, period=kper) + else + varid = nc_vars%varid(idt%mf6varname) + end if call load_integer1d_type(int1d, mf6_input, mshape, idt, nc_vars, & varid, input_fname) end if @@ -100,7 +105,7 @@ subroutine nc_array_load_int2d(int2d, mshape, idt, mf6_input, nc_vars, & call load_integer2d_layered(int2d, mf6_input, mshape, idt, nc_vars, & input_fname) else - varid = nc_vars%varid(idt%tagname) + varid = nc_vars%varid(idt%mf6varname) call load_integer2d_type(int2d, mf6_input, mshape, idt, nc_vars, & varid, input_fname) end if @@ -130,7 +135,7 @@ subroutine nc_array_load_int3d(int3d, mshape, idt, mf6_input, nc_vars, & call load_integer3d_layered(int3d, mf6_input, mshape, idt, nc_vars, & input_fname) else - varid = nc_vars%varid(idt%tagname) + varid = nc_vars%varid(idt%mf6varname) call load_integer3d_type(int3d, mf6_input, mshape, idt, nc_vars, & varid, input_fname) end if @@ -142,7 +147,7 @@ end subroutine nc_array_load_int3d !> @brief Load NetCDF double 1D array !< subroutine nc_array_load_dbl1d(dbl1d, mshape, idt, mf6_input, nc_vars, & - input_fname, iout) + input_fname, iout, kper, iaux) real(DP), dimension(:), pointer, contiguous, intent(in) :: dbl1d integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape !< model shape type(InputParamDefinitionType), intent(in) :: idt !< input data type object describing this record @@ -150,19 +155,35 @@ subroutine nc_array_load_dbl1d(dbl1d, mshape, idt, mf6_input, nc_vars, & type(NCPackageVarsType), pointer, intent(in) :: nc_vars character(len=*), intent(in) :: input_fname integer(I4B), intent(in) :: iout + integer(I4B), optional, intent(in) :: kper + integer(I4B), optional, intent(in) :: iaux ! -- local integer(I4B) :: varid logical(LGP) :: layered ! - layered = (idt%layered .and. is_layered(nc_vars%grid)) + if (present(kper)) then + layered = (kper > 0 .and. is_layered(nc_vars%grid)) + else + layered = (idt%layered .and. is_layered(nc_vars%grid)) + end if ! if (layered) then - call load_double1d_layered(dbl1d, mf6_input, mshape, idt, nc_vars, & - input_fname) + if (present(kper)) then + call load_double1d_layered_spd(dbl1d, mf6_input, mshape, idt, nc_vars, & + kper, input_fname, iaux) + else + call load_double1d_layered(dbl1d, mf6_input, mshape, idt, nc_vars, & + input_fname) + end if else - varid = nc_vars%varid(idt%tagname) - call load_double1d_type(dbl1d, mf6_input, mshape, idt, nc_vars, & - varid, input_fname) + if (present(kper)) then + call load_double1d_spd(dbl1d, mf6_input, mshape, idt, nc_vars, & + kper, input_fname, iaux) + else + varid = nc_vars%varid(idt%mf6varname) + call load_double1d_type(dbl1d, mf6_input, mshape, idt, nc_vars, & + varid, input_fname) + end if end if ! ! -- return @@ -190,7 +211,7 @@ subroutine nc_array_load_dbl2d(dbl2d, mshape, idt, mf6_input, nc_vars, & call load_double2d_layered(dbl2d, mf6_input, mshape, idt, nc_vars, & input_fname) else - varid = nc_vars%varid(idt%tagname) + varid = nc_vars%varid(idt%mf6varname) call load_double2d_type(dbl2d, mf6_input, mshape, idt, nc_vars, & varid, input_fname) end if @@ -220,7 +241,7 @@ subroutine nc_array_load_dbl3d(dbl3d, mshape, idt, mf6_input, nc_vars, & call load_double3d_layered(dbl3d, mf6_input, mshape, idt, nc_vars, & input_fname) else - varid = nc_vars%varid(idt%tagname) + varid = nc_vars%varid(idt%mf6varname) call load_double3d_type(dbl3d, mf6_input, mshape, idt, nc_vars, & varid, input_fname) end if @@ -305,7 +326,7 @@ subroutine load_integer1d_layered(int1d, mf6_input, mshape, idt, nc_vars, & ncpl = product(layer_shape) index_start = 1 do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) index_stop = index_start + ncpl - 1 int1d_ptr(1:ncpl) => int1d(index_start:index_stop) call nf_verify(nf90_get_var(nc_vars%ncid, varid, int1d_ptr), & @@ -379,7 +400,7 @@ subroutine load_integer2d_layered(int2d, mf6_input, mshape, idt, nc_vars, & call get_layered_shape(mshape, nlay, layer_shape) ncpl = layer_shape(1) do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) int1d_ptr(1:ncpl) => int2d(1:ncpl, k) call nf_verify(nf90_get_var(nc_vars%ncid, varid, int1d_ptr), & nc_vars%nc_fname) @@ -435,7 +456,7 @@ subroutine load_integer3d_layered(int3d, mf6_input, mshape, idt, nc_vars, & ncpl = product(layer_shape) ! do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) index_stop = index_start + ncpl - 1 int1d_ptr(1:ncpl) => int3d(:, :, k:k) call nf_verify(nf90_get_var(nc_vars%ncid, varid, int1d_ptr), & @@ -498,6 +519,78 @@ subroutine load_double1d_type(dbl1d, mf6_input, mshape, idt, nc_vars, & return end subroutine load_double1d_type + !> @brief load type 1d double + !< + subroutine load_double1d_spd(dbl1d, mf6_input, mshape, idt, nc_vars, & + iper, input_fname, iaux) + use ConstantsModule, only: DNODATA + ! -- dummy + real(DP), dimension(:), contiguous, pointer, intent(in) :: dbl1d + type(ModflowInputType), intent(in) :: mf6_input + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape + type(InputParamDefinitionType), intent(in) :: idt + type(NCPackageVarsType), pointer, intent(in) :: nc_vars + integer(I4B), intent(in) :: iper + character(len=*), intent(in) :: input_fname + integer(I4B), optional, intent(in) :: iaux + ! -- local + real(DP), dimension(:, :, :), contiguous, pointer :: dbl3d + integer(I4B) :: nvals, varid + integer(I4B) :: n, i, j, k + ! + ! -- initialize + nvals = 0 + ! + ! -- set varid + if (present(iaux)) then + varid = nc_vars%varid(idt%mf6varname, period=iper, iaux=iaux) + else + varid = nc_vars%varid(idt%mf6varname, period=iper) + end if + ! + if (idt%shape == 'NODES') then + ! TODO future support + write (errmsg, '(a)') & + 'IDM NetCDF load_double1d_spd NODES var shape not supported => '// & + trim(idt%tagname) + call store_error(errmsg) + call store_error_filename(input_fname) + else if (idt%shape == 'NCPL' .or. idt%shape == 'NAUX NCPL') then + + if (size(mshape) == 3) then + ! + allocate (dbl3d(mshape(3), mshape(2), mshape(1))) + ! + call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl3d), & + nc_vars%nc_fname) + ! + n = 0 + do k = 1, size(dbl3d, dim=3) + do i = 1, size(dbl3d, dim=2) + do j = 1, size(dbl3d, dim=1) + if (n < size(dbl1d)) then + n = n + 1 + else + n = 1 + end if + if (dbl3d(j, i, k) /= DNODATA) then + dbl1d(n) = dbl3d(j, i, k) + end if + end do + end do + end do + + else if (size(mshape) == 2) then + ! TODO + write (errmsg, '(a)') & + 'IDM NetCDF load_double1d_spd DISV model not supported => '// & + trim(idt%tagname) + call store_error(errmsg) + call store_error_filename(input_fname) + end if + end if + end subroutine load_double1d_spd + !> @brief load type 1d double layered !< subroutine load_double1d_layered(dbl1d, mf6_input, mshape, idt, nc_vars, & @@ -524,7 +617,7 @@ subroutine load_double1d_layered(dbl1d, mf6_input, mshape, idt, nc_vars, & ncpl = product(layer_shape) ! do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) index_stop = index_start + ncpl - 1 dbl1d_ptr(1:ncpl) => dbl1d(index_start:index_stop) call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl1d_ptr), & @@ -536,6 +629,56 @@ subroutine load_double1d_layered(dbl1d, mf6_input, mshape, idt, nc_vars, & return end subroutine load_double1d_layered + !> @brief load type 1d double layered + !< + subroutine load_double1d_layered_spd(dbl1d, mf6_input, mshape, idt, nc_vars, & + iper, input_fname, iaux) + use ConstantsModule, only: DNODATA + ! -- dummy + real(DP), dimension(:), contiguous, pointer, intent(in) :: dbl1d + type(ModflowInputType), intent(in) :: mf6_input + integer(I4B), dimension(:), contiguous, pointer, intent(in) :: mshape + type(InputParamDefinitionType), intent(in) :: idt + type(NCPackageVarsType), pointer, intent(in) :: nc_vars + integer(I4B), intent(in) :: iper + character(len=*), intent(in) :: input_fname + integer(I4B), optional, intent(in) :: iaux + ! -- local + integer(I4B), dimension(:), allocatable :: layer_shape + integer(I4B) :: nlay, varid + integer(I4B) :: k, n, ncpl + integer(I4B) :: index_start, index_stop + real(DP), dimension(:), contiguous, pointer :: dbl1d_ptr + ! + call get_layered_shape(mshape, nlay, layer_shape) + ! + ncpl = product(layer_shape) + allocate (dbl1d_ptr(ncpl)) + ! + do k = 1, nlay + index_start = 1 + index_stop = index_start + ncpl - 1 + if (present(iaux)) then + varid = nc_vars%varid(idt%mf6varname, layer=k, period=iper, iaux=iaux) + else + varid = nc_vars%varid(idt%mf6varname, layer=k, period=iper) + end if + call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl1d_ptr), & + nc_vars%nc_fname) + do n = 1, ncpl + if (dbl1d_ptr(n) /= DNODATA) then + dbl1d(n) = dbl1d_ptr(n) + end if + end do + end do + ! + ! -- cleanup + deallocate (dbl1d_ptr) + ! + ! -- return + return + end subroutine load_double1d_layered_spd + !> @brief load type 2d double !< subroutine load_double2d_type(dbl2d, mf6_input, mshape, idt, nc_vars, varid, & @@ -598,7 +741,7 @@ subroutine load_double2d_layered(dbl2d, mf6_input, mshape, idt, nc_vars, & call get_layered_shape(mshape, nlay, layer_shape) ncpl = layer_shape(1) do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) dbl1d_ptr(1:ncpl) => dbl2d(1:ncpl, k) call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl1d_ptr), & nc_vars%nc_fname) @@ -653,7 +796,7 @@ subroutine load_double3d_layered(dbl3d, mf6_input, mshape, idt, nc_vars, & ncpl = product(layer_shape) index_start = 1 do k = 1, nlay - varid = nc_vars%varid(idt%tagname, k) + varid = nc_vars%varid(idt%mf6varname, layer=k) index_stop = index_start + ncpl - 1 dbl1d_ptr(1:ncpl) => dbl3d(:, :, k:k) call nf_verify(nf90_get_var(nc_vars%ncid, varid, dbl1d_ptr), & diff --git a/src/Utilities/Idm/netcdf/NCContextBuild.f90 b/src/Utilities/Idm/netcdf/NCContextBuild.f90 index 12ebf0d0f3a..7ded2b6b591 100644 --- a/src/Utilities/Idm/netcdf/NCContextBuild.f90 +++ b/src/Utilities/Idm/netcdf/NCContextBuild.f90 @@ -73,11 +73,13 @@ subroutine add_package_var(modeltype, modelname, nc_vars, input_name, varid, & character(len=NETCDF_ATTR_STRLEN) :: input_str character(len=LENCOMPONENTNAME) :: c_name, sc_name character(len=LINELENGTH) :: mempath, varname - integer(I4B) :: layer, mf6_layer + integer(I4B) :: layer, period, iaux, mf6_layer, mf6_period, mf6_iaux logical(LGP) :: success ! ! -- initialize layer = -1 + period = -1 + iaux = -1 varname = '' c_name = '' sc_name = '' @@ -104,8 +106,20 @@ subroutine add_package_var(modeltype, modelname, nc_vars, input_name, varid, & layer = mf6_layer end if ! + ! -- check for optional period attribute + if (nf90_get_att(nc_vars%ncid, varid, & + 'modflow6_iper', mf6_period) == NF90_NOERR) then + period = mf6_period + end if + ! + ! -- check for optional period attribute + if (nf90_get_att(nc_vars%ncid, varid, & + 'modflow6_iaux', mf6_iaux) == NF90_NOERR) then + iaux = mf6_iaux + end if + ! ! -- add the variable to netcdf description - call nc_vars%add(sc_name, varname, layer, varid) + call nc_vars%add(sc_name, varname, layer, period, iaux, varid) else errmsg = 'NetCDF variable invalid modflow6_input attribute: "'// & trim(input_str)//'".' diff --git a/src/Utilities/Idm/netcdf/NCFileVars.f90 b/src/Utilities/Idm/netcdf/NCFileVars.f90 index ccecd844d3f..d7cd82cbcc9 100644 --- a/src/Utilities/Idm/netcdf/NCFileVars.f90 +++ b/src/Utilities/Idm/netcdf/NCFileVars.f90 @@ -21,9 +21,7 @@ module NCFileVarsModule !< type :: NCPackageVarsType character(len=LENMODELNAME) :: modelname !< name of model - character(len=LINELENGTH), dimension(:), allocatable :: tagnames !< variable tag name - integer(I4B), dimension(:), allocatable :: layers !< variable layer - integer(I4B), dimension(:), allocatable :: varids !< netcdf file variable id + type(ListType) :: nc_vars character(len=LINELENGTH), pointer :: grid => null() !< grid type character(len=LINELENGTH), pointer :: nc_fname => null() !< netcdf filename integer(I4B), pointer :: ncid => null() !< netcdf file handle @@ -39,6 +37,8 @@ module NCFileVarsModule character(LINELENGTH) :: pkgname !< package name character(LINELENGTH) :: tagname !< tag name integer(I4B) :: layer !< variable layer + integer(I4B) :: period !< variable period + integer(I4B) :: iaux !< variable aux index integer(I4B) :: varid !< NC file variable id contains end type NCFileMf6VarType @@ -47,15 +47,12 @@ module NCFileVarsModule !< type :: NCFileVarsType type(ListType) :: mf6invar !< list of modflow 6 input variables in netcdf file - character(len=LINELENGTH), dimension(:), allocatable :: pkgnames !< packages in file - integer(I4B), dimension(:), allocatable :: pkgcounts !< variable counts in each package character(len=LINELENGTH), pointer :: grid => null() !< grid type character(len=LINELENGTH), pointer :: nc_fname => null() !< netcdf filename integer(I4B), pointer :: ncid => null() !< netcdf file handle contains procedure :: init => fv_init procedure :: add => fv_add - procedure :: get => fv_get procedure :: destroy => fv_destroy procedure :: create_varlists end type NCFileVarsType @@ -74,56 +71,71 @@ subroutine ncvars_init(this, modelname) ! -- set modelname this%modelname = modelname ! - ! -- allocate empty arrays - allocate (this%tagnames(0)) - allocate (this%layers(0)) - allocate (this%varids(0)) - ! ! -- return return end subroutine ncvars_init !> @brief return a netcdf variable id for a package tagname !< - function ncvars_varid(this, tagname, layer) result(varid) + function ncvars_varid(this, tagname, layer, period, iaux) result(varid) ! -- modules ! -- dummy class(NCPackageVarsType) :: this character(len=*), intent(in) :: tagname integer(I4B), optional :: layer + integer(I4B), optional :: period + integer(I4B), optional :: iaux ! -- return integer(I4B) :: varid ! -- local - integer(I4B) :: n, l + integer(I4B) :: n, l, p, a + class(NCFileMf6VarType), pointer :: nc_var ! ! -- initialize varid = -1 l = -1 + p = -1 + a = -1 ! ! -- set search layer if provided if (present(layer)) then l = layer end if ! - do n = 1, size(this%tagnames) - if (this%tagnames(n) == tagname .and. & - this%layers(n) == l) then - varid = this%varids(n) - exit + ! -- set search period if provided + if (present(period)) then + p = period + end if + ! -- set search iaux if provided + if (present(iaux)) then + a = iaux + end if + ! + do n = 1, this%nc_vars%Count() + nc_var => ncvar_get(this%nc_vars, n) + if (nc_var%tagname == tagname .and. & + nc_var%layer == l .and. & + nc_var%period == p .and. & + nc_var%iaux == a) then + varid = nc_var%varid end if end do ! ! -- set error and exit if variable not in NetCDF input if (varid == -1) then if (this%nc_fname /= '') then - if (l == -1) then - write (errmsg, '(a)') & - 'NetCDF variable not found, tagname="'//trim(tagname)//'".' - else - write (errmsg, '(a,i0,a)') & - 'NetCDF variable not found, tagname="'//trim(tagname)// & - '", layer=', l, '.' + write (errmsg, '(a)') & + 'NetCDF variable not found, tagname="'//trim(tagname)//'"' + if (present(layer)) then + write (errmsg, '(a,i0)') trim(errmsg)//', ilayer=', layer end if + if (present(period)) then + write (errmsg, '(a,i0)') trim(errmsg)//', period=', period + end if + if (present(iaux)) then + write (errmsg, '(a,i0)') trim(errmsg)//', iaux=', iaux + end if + write (errmsg, '(a)') trim(errmsg)//'.' call store_error(errmsg) call store_error_filename(this%nc_fname) else @@ -145,11 +157,17 @@ subroutine ncvars_destroy(this) ! -- dummy class(NCPackageVarsType) :: this ! -- local + class(NCFileMf6VarType), pointer :: nc_var + integer(I4B) :: n + ! + ! -- deallocate allocated memory + do n = 1, this%nc_vars%Count() + nc_var => ncvar_get(this%nc_vars, n) + deallocate (nc_var) + nullify (nc_var) + end do ! - ! -- deallocate dynamic arrays - if (allocated(this%tagnames)) deallocate (this%tagnames) - if (allocated(this%layers)) deallocate (this%layers) - if (allocated(this%varids)) deallocate (this%varids) + call this%nc_vars%Clear() ! ! -- return return @@ -184,10 +202,6 @@ subroutine fv_init(this, modelname, nc_fname, ncid, grid) call mem_allocate(this%nc_fname, ilen, 'NETCDF_FNAME', mempath) call mem_allocate(this%ncid, 'NCID', mempath) ! - ! -- allocate local memory - allocate (this%pkgnames(0)) - allocate (this%pkgcounts(0)) - ! ! -- set this%grid = grid this%nc_fname = nc_fname @@ -199,7 +213,7 @@ end subroutine fv_init !> @brief add netcdf modflow6 input variable to list !< - subroutine fv_add(this, pkgname, tagname, layer, varid) + subroutine fv_add(this, pkgname, tagname, layer, period, iaux, varid) ! -- modules use ArrayHandlersModule, only: expandarray ! -- dummy variables @@ -207,39 +221,20 @@ subroutine fv_add(this, pkgname, tagname, layer, varid) character(len=*), intent(in) :: pkgname character(len=*), intent(in) :: tagname integer(I4B), intent(in) :: layer + integer(I4B), intent(in) :: period + integer(I4B), intent(in) :: iaux integer(I4B), intent(in) :: varid ! -- local variables class(NCFileMf6VarType), pointer :: invar class(*), pointer :: obj - integer(I4B) :: n, pidx - ! - ! -- initialize package index - pidx = 0 - ! - do n = 1, size(this%pkgnames) - if (this%pkgnames(n) == pkgname) then - pidx = n - exit - end if - end do - ! - ! -- add new package type - if (pidx == 0) then - call expandarray(this%pkgnames) - call expandarray(this%pkgcounts) - pidx = size(this%pkgnames) - this%pkgnames(pidx) = pkgname - this%pkgcounts(pidx) = 0 - end if - ! - ! --increment pkgcount - this%pkgcounts(pidx) = this%pkgcounts(pidx) + 1 ! ! -- add mf6 variable to file list allocate (invar) invar%pkgname = pkgname invar%tagname = tagname invar%layer = layer + invar%period = period + invar%iaux = iaux invar%varid = varid ! obj => invar @@ -249,32 +244,6 @@ subroutine fv_add(this, pkgname, tagname, layer, varid) return end subroutine fv_add - !> @brief get modflow6 input variable description at position idx - !< - function fv_get(this, idx) result(res) - ! -- dummy variables - class(NCFileVarsType) :: this - integer(I4B), intent(in) :: idx !< package number - class(NCFileMf6VarType), pointer :: res - ! -- local variables - class(*), pointer :: obj - ! - ! -- initialize res - res => null() - ! - ! -- get the package from the list - obj => this%mf6invar%GetItem(idx) - if (associated(obj)) then - select type (obj) - class is (NCFileMf6VarType) - res => obj - end select - end if - ! - ! -- return - return - end function fv_get - !> @brief destroy netcdf model variable description type !< subroutine fv_destroy(this) @@ -285,16 +254,13 @@ subroutine fv_destroy(this) integer(I4B) :: n ! do n = 1, this%mf6invar%Count() - invar => this%get(n) + invar => ncvar_get(this%mf6invar, n) deallocate (invar) nullify (invar) end do ! call this%mf6invar%Clear() ! - if (allocated(this%pkgnames)) deallocate (this%pkgnames) - if (allocated(this%pkgcounts)) deallocate (this%pkgcounts) - ! return end subroutine fv_destroy @@ -307,48 +273,27 @@ subroutine create_varlists(this, modelname, pkgname, nc_vars) character(len=*), intent(in) :: modelname character(len=*), intent(in) :: pkgname type(NCPackageVarsType), pointer, intent(inout) :: nc_vars - integer(I4B) :: n, cnt, pidx + integer(I4B) :: n ! -- local - class(NCFileMf6VarType), pointer :: invar - ! - ! -- initialize - cnt = 0 - pidx = 0 - ! - ! -- deallocate incoming lists - if (allocated(nc_vars%tagnames)) deallocate (nc_vars%tagnames) - if (allocated(nc_vars%layers)) deallocate (nc_vars%layers) - if (allocated(nc_vars%varids)) deallocate (nc_vars%varids) + class(NCFileMf6VarType), pointer :: invar, nc_var + class(*), pointer :: obj ! - do n = 1, size(this%pkgnames) - if (this%pkgnames(n) == pkgname) then - pidx = n - exit + do n = 1, this%mf6invar%count() + invar => ncvar_get(this%mf6invar, n) + if (invar%pkgname == pkgname) then + ! -- create package variable description + allocate (nc_var) + nc_var%pkgname = invar%pkgname + nc_var%tagname = invar%tagname + nc_var%layer = invar%layer + nc_var%period = invar%period + nc_var%iaux = invar%iaux + nc_var%varid = invar%varid + obj => nc_var + call nc_vars%nc_vars%Add(obj) end if end do ! - if (pidx > 0) then - ! -- package has NCFile variables - ! - allocate (nc_vars%tagnames(this%pkgcounts(pidx))) - allocate (nc_vars%layers(this%pkgcounts(pidx))) - allocate (nc_vars%varids(this%pkgcounts(pidx))) - ! - do n = 1, this%mf6invar%count() - invar => this%get(n) - if (invar%pkgname == pkgname) then - cnt = cnt + 1 - nc_vars%tagnames(cnt) = invar%tagname - nc_vars%layers(cnt) = invar%layer - nc_vars%varids(cnt) = invar%varid - end if - end do - else - allocate (nc_vars%tagnames(0)) - allocate (nc_vars%layers(0)) - allocate (nc_vars%varids(0)) - end if - ! ! -- set modelname nc_vars%modelname = modelname ! @@ -361,4 +306,30 @@ subroutine create_varlists(this, modelname, pkgname, nc_vars) return end subroutine create_varlists + !> @brief get modflow6 input variable description at position idx + !< + function ncvar_get(nc_vars, idx) result(res) + ! -- dummy variables + type(ListType) :: nc_vars + integer(I4B), intent(in) :: idx !< package number + class(NCFileMf6VarType), pointer :: res + ! -- local variables + class(*), pointer :: obj + ! + ! -- initialize res + res => null() + ! + ! -- get the package from the list + obj => nc_vars%GetItem(idx) + if (associated(obj)) then + select type (obj) + class is (NCFileMf6VarType) + res => obj + end select + end if + ! + ! -- return + return + end function ncvar_get + end module NCFileVarsModule diff --git a/src/mf6core.f90 b/src/mf6core.f90 index 44fd1a14a2b..1c3cedbf8e7 100644 --- a/src/mf6core.f90 +++ b/src/mf6core.f90 @@ -482,6 +482,7 @@ subroutine Mf6PrepareTimestep() use SimModule, only: converge_reset use SimVariablesModule, only: isim_mode use IdmLoadModule, only: idm_rp + use SourceLoadModule, only: export_post_prepare ! -- local variables class(BaseModelType), pointer :: mp => null() class(BaseExchangeType), pointer :: ep => null() @@ -567,6 +568,9 @@ subroutine Mf6PrepareTimestep() call sp%sln_dt() end do ! + ! -- update exports + call export_post_prepare() + ! ! -- set time step call tdis_set_timestep()