diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index ed25a84c..a2af53e8 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -175,10 +175,11 @@ jobs: pip list - name: Run available domains with PRMS and convert csv output to NetCDF - working-directory: test_data/scripts + working-directory: test_data/generate run: | - pytest -v -n=auto --durations=0 test_run_domains.py - pytest -v -n=auto --durations=0 test_nc_domains.py + pytest -v -n=auto --durations=0 run_prms_domains.py + pytest -v -n=auto --durations=0 convert_prms_output_to_nc.py + pytest -v -n=auto --durations=0 remove_prms_csvs.py - name: List all NetCDF files in test_data directory working-directory: test_data @@ -198,13 +199,13 @@ jobs: - name: Upload test results if: always() - uses: actions/upload-artifact@v2 + uses: actions/upload-artifact@v3 with: name: Test results for ${{ runner.os }}-${{ matrix.python-version }} path: ./autotest/pytest.xml - name: Upload code coverage to Codecov - uses: codecov/codecov-action@v2.1.0 + uses: codecov/codecov-action@v3 with: file: ./autotest/coverage.xml # flags: unittests diff --git a/.gitignore b/.gitignore index 8f0e88e1..4c7f1c3f 100644 --- a/.gitignore +++ b/.gitignore @@ -61,12 +61,15 @@ generated/ # example output data and gis data examples/0*/ +examples/snow_errors +examples/runoff_errors examples/pynhm_gis examples/pywatershed_gis examples/model_loop_custom_output pywatershed/data/pynhm_gis pywatershed/data/pywatershed_gis + # graphics *.png *.svg diff --git a/autotest/test_control.py b/autotest/test_control.py index 4b215141..946425c1 100644 --- a/autotest/test_control.py +++ b/autotest/test_control.py @@ -80,7 +80,7 @@ def test_control_simple(control_simple): year = year if month >= 10 else year - 1 wy_start = np.datetime64(f"{year}-10-01") dowy = (current_time - wy_start).astype("timedelta64[D]") - assert dowy == control_simple.current_dowy + assert dowy == (control_simple.current_dowy - 1) prev_time = control_simple.current_time diff --git a/autotest/test_model.py b/autotest/test_model.py index 0f15ad56..8cc1bdd9 100644 --- a/autotest/test_model.py +++ b/autotest/test_model.py @@ -287,18 +287,18 @@ def test_model(domain, model_args, tmp_path): 9: { "PRMSChannel": { "seg_outflow": { - "drb_2yr": 1430.6364027142613, - "hru_1": 13.416914151483681, - "ucb_2yr": 1694.5412856707849, + "drb_2yr": 1517.232887980279, + "hru_1": 13.696918669514927, + "ucb_2yr": 1694.5697712423928, }, }, }, 99: { "PRMSChannel": { "seg_outflow": { - "drb_2yr": 1588.1444684289775, - "hru_1": 19.596412903692578, - "ucb_2yr": 407.2200022510677, + "drb_2yr": 2350.499659332901, + "hru_1": 22.874414994530095, + "ucb_2yr": 733.2293013532435, }, }, }, diff --git a/autotest/test_nhm_self_drive.py b/autotest/test_nhm_self_drive.py index c0a2ecc5..f8b734ba 100644 --- a/autotest/test_nhm_self_drive.py +++ b/autotest/test_nhm_self_drive.py @@ -22,10 +22,14 @@ def test_drive_indiv_process(domain, tmp_path): - """Use output from a full NHM run to drive each of the indiv processes - separately: self-driving + """Output of a full pywatershed NHM drives indiv process models separately + + The results from the full model should be consistent with the results from + the individual models, else there is likely something wrong with the + full model. """ - # Full NHM output + + # Run a full pws NHM to use its output to drive individual processes nhm_output_dir = pl.Path(tmp_path) / "nhm_output" params = pws.parameters.PrmsParameters.load(domain["param_file"]) @@ -44,10 +48,9 @@ def test_drive_indiv_process(domain, tmp_path): nhm.run(finalize=True) del nhm, params, control - # individual process models + # run individual process models for proc in nhm_processes: - # proc = pws.PRMSRunoff # TODO: fix this one ASAP - if proc in [pws.PRMSSolarGeometry, pws.PRMSAtmosphere, pws.PRMSRunoff]: + if proc in [pws.PRMSSolarGeometry, pws.PRMSAtmosphere]: # These are not driven by outputs of above, only external outputs # or known/static inputs continue @@ -82,10 +85,11 @@ def test_drive_indiv_process(domain, tmp_path): ans = xr.open_dataset(nhm_output_dir / f"{vv}.nc") # Leaving the commented to diagnose what PRMSRunoff later. - # try: - xr.testing.assert_allclose(res, ans) - # except: - # print(vv) + try: + xr.testing.assert_allclose(res, ans) + except: + print(vv, abs(res - ans).max()) + print(vv, (abs(res - ans) / ans).max()) del res, ans diff --git a/autotest/test_preprocess_csv.py b/autotest/test_preprocess_csv.py index 0364411c..8ce8bfb1 100644 --- a/autotest/test_preprocess_csv.py +++ b/autotest/test_preprocess_csv.py @@ -7,6 +7,8 @@ from pywatershed import CsvFile +# these CSV files are protected from deletion in CI by +# test_data/scripts/test_remove_csvs.py csv_test_vars = ["hru_ppt", "intcp_stor", "potet", "gwres_stor"] diff --git a/autotest/test_prms_channel.py b/autotest/test_prms_channel.py index 299db4a7..5160e816 100644 --- a/autotest/test_prms_channel.py +++ b/autotest/test_prms_channel.py @@ -2,11 +2,16 @@ import pytest +from pywatershed.base.adapter import adapter_factory from pywatershed.base.control import Control from pywatershed.base.parameters import Parameters from pywatershed.hydrology.prms_channel import PRMSChannel, has_prmschannel_f from pywatershed.parameters import PrmsParameters -from pywatershed.utils.netcdf_utils import NetCdfCompare +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? +compare_output_files = False +rtol = atol = 1.0e-7 fail_fast = False @@ -67,52 +72,39 @@ def test_compare_prms( budget_type="error", calc_method=calc_method, ) - nc_parent = tmp_path / domain["domain_name"] - channel.initialize_netcdf(nc_parent) - # test that init netcdf twice raises a warning - with pytest.warns(UserWarning): + + if compare_output_files: + nc_parent = tmp_path / domain["domain_name"] channel.initialize_netcdf(nc_parent) + # test that init netcdf twice raises a warning + with pytest.warns(UserWarning): + channel.initialize_netcdf(nc_parent) + + else: + answers = {} + for var in PRMSChannel.get_variables(): + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) for istep in range(control.n_times): control.advance() - channel.advance() - channel.calculate(float(istep)) - channel.output() + if not compare_output_files: + compare_in_memory(channel, answers, atol=atol, rtol=rtol) channel.finalize() - output_compare = {} - - for key in PRMSChannel.get_variables(): - base_nc_path = output_dir / f"{key}.nc" - compare_nc_path = tmp_path / domain["domain_name"] / f"{key}.nc" - # PRMS does not output the storage change in the channel - if not base_nc_path.exists(): - continue - output_compare[key] = (base_nc_path, compare_nc_path) - - assert_error = False - for key, (base, compare) in output_compare.items(): - print(f"\nbase_nc_path: {base}") - print(f"compare_nc_path: {compare}") - success, diff = NetCdfCompare(base, compare).compare() - if not success: - print( - f"comparison for {key} failed: " - + f"maximum error {diff[key][0]} " - + f"(maximum allowed error {diff[key][1]}) " - + f"in column {diff[key][2]}" - ) - assert_error = True - if fail_fast: - assert False - - else: - print(f"comparison for {key} passed") - - assert not assert_error, "comparison failed" + if compare_output_files: + compare_netcdfs( + PRMSChannel.get_variables(), + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return diff --git a/autotest/test_prms_groundwater.py b/autotest/test_prms_groundwater.py index a3381fe0..f0209d8d 100644 --- a/autotest/test_prms_groundwater.py +++ b/autotest/test_prms_groundwater.py @@ -2,14 +2,15 @@ import pytest -from pywatershed.base.control import Control -from pywatershed.base.parameters import Parameters -from pywatershed.hydrology.prms_groundwater import ( - PRMSGroundwater, - has_prmsgroundwater_f, -) +from pywatershed import Control, Parameters, PRMSGroundwater +from pywatershed.base.adapter import adapter_factory +from pywatershed.hydrology.prms_groundwater import has_prmsgroundwater_f from pywatershed.parameters import PrmsParameters -from pywatershed.utils.netcdf_utils import NetCdfCompare +from utils_compare import compare_in_memory, compare_netcdfs + +# compare in memory (faster) or full output files? +compare_output_files = False +rtol = atol = 1.0e-13 calc_methods = ("numpy", "numba", "fortran") params = ("params_sep", "params_one") @@ -48,7 +49,6 @@ def test_compare_prms( tmp_path = pl.Path(tmp_path) - # load csv files into dataframes output_dir = domain["prms_output_dir"] input_variables = {} for key in PRMSGroundwater.get_inputs(): @@ -63,49 +63,36 @@ def test_compare_prms( budget_type="error", calc_method=calc_method, ) - nc_parent = tmp_path / domain["domain_name"] - gw.initialize_netcdf(nc_parent) - - output_compare = {} - vars_compare = ( - "gwres_flow", - "gwres_sink", - "gwres_stor", - "ssr_to_gw", - "soil_to_gw", - ) - for key in PRMSGroundwater.get_variables(): - if key not in vars_compare: - continue - base_nc_path = output_dir / f"{key}.nc" - compare_nc_path = tmp_path / domain["domain_name"] / f"{key}.nc" - output_compare[key] = (base_nc_path, compare_nc_path) - print(f"base_nc_path: {base_nc_path}") - print(f"compare_nc_path: {compare_nc_path}") + if compare_output_files: + nc_parent = tmp_path / domain["domain_name"] + gw.initialize_netcdf(nc_parent) + else: + answers = {} + for var in PRMSGroundwater.get_variables(): + var_pth = output_dir / f"{var}.nc" + answers[var] = adapter_factory( + var_pth, variable_name=var, control=control + ) for istep in range(control.n_times): control.advance() - gw.advance() - gw.calculate(float(istep)) - gw.output() + if not compare_output_files: + compare_in_memory(gw, answers, atol=atol, rtol=rtol) + gw.finalize() - assert_error = False - for key, (base, compare) in output_compare.items(): - success, diff = NetCdfCompare(base, compare).compare() - if not success: - print( - f"comparison for {key} failed: " - + f"maximum error {diff[key][0]} " - + f"(maximum allowed error {diff[key][1]}) " - + f"in column {diff[key][2]}" - ) - assert_error = True - assert not assert_error, "comparison failed" + if compare_output_files: + compare_netcdfs( + PRMSGroundwater.get_variables(), + tmp_path / domain["domain_name"], + output_dir, + atol=atol, + rtol=rtol, + ) return diff --git a/autotest/test_prms_runoff.py b/autotest/test_prms_runoff.py index 88a25517..f0fe0de6 100644 --- a/autotest/test_prms_runoff.py +++ b/autotest/test_prms_runoff.py @@ -51,6 +51,15 @@ def test_compare_prms( "sroff", "dprst_evap_hru", "hru_impervevap", + "dprst_insroff_hru", + "dprst_stor_hru", + "contrib_fraction", + "hru_sroffp", + "hru_sroffi", + # "hru_impervstor_change", + "dprst_sroff_hru", + "dprst_insroff_hru", + # "dprst_stor_hru_change", ] output_dir = domain["prms_output_dir"] @@ -72,7 +81,7 @@ def test_compare_prms( parameters=parameters, calc_method=calc_method, **input_variables, - budget_type="warn", # intermittent errors currently + budget_type="error", ) all_success = True @@ -90,7 +99,7 @@ def test_compare_prms( failfast = True detailed = True if check: - atol = 1.0e-5 + atol = 1.0e-10 success = check_timestep_results( runoff, istep, ans, atol, detailed ) @@ -113,7 +122,7 @@ def check_timestep_results(storageunit, istep, ans, atol, detailed=False): for key in ans.keys(): a1 = ans[key].current a2 = storageunit[key] - success = np.isclose(a1, a2, atol=atol).all() + success = np.isclose(a1, a2, atol=atol, rtol=atol).all() if not success: all_success = False diff = a1 - a2 diff --git a/autotest/test_prms_solar_geom.py b/autotest/test_prms_solar_geom.py index 4d40cd9a..f2a4560c 100644 --- a/autotest/test_prms_solar_geom.py +++ b/autotest/test_prms_solar_geom.py @@ -31,7 +31,6 @@ def parameters(domain, request): return params -@pytest.mark.xfail @pytest.mark.parametrize( "from_prms_file", (True, False), ids=("from_prms_file", "compute") ) diff --git a/autotest/utils_compare.py b/autotest/utils_compare.py new file mode 100644 index 00000000..dd379184 --- /dev/null +++ b/autotest/utils_compare.py @@ -0,0 +1,121 @@ +import pathlib as pl + +import numpy as np +import xarray as xr + +import pywatershed as pws + + +def assert_allclose( + actual: np.ndarray, + desired: np.ndarray, + rtol: float = 1.0e-15, + atol: float = 1.0e-15, + equal_nan: bool = True, + strict: bool = False, + also_check_w_np: bool = True, + error_message: str = "Comparison unsuccessful (default message)", +): + """Reinvent np.testing.assert_allclose to get useful diagnostincs in debug + + Args: + actual: Array obtained. + desired: Array desired. + rtol: Relative tolerance. + atol: Absolute tolerance. + equal_nan: If True, NaNs will compare equal. + strict: If True, raise an ``AssertionError`` when either the shape or + the data type of the arguments does not match. The special + handling of scalars mentioned in the Notes section is disabled. + also_check_w_np: first check using np.testing.assert_allclose using + the same options. + """ + + if also_check_w_np: + np.testing.assert_allclose( + actual, + desired, + rtol=rtol, + atol=atol, + equal_nan=equal_nan, + # strict=strict, # to add to newer versions of numpy + ) + + if strict: + assert actual.shape == desired.shape + assert isinstance(actual, type(desired)) + assert isinstance(desired, type(actual)) + + if equal_nan: + actual_nan = np.where(np.isnan(actual), True, False) + desired_nan = np.where(np.isnan(desired), True, False) + assert (actual_nan == desired_nan).all() + + abs_diff = abs(actual - desired) + rel_abs_diff = abs_diff / desired + + abs_close = abs_diff < atol + rel_close = rel_abs_diff < rtol + rel_close = np.where(np.isnan(rel_close), False, rel_close) + + close = abs_close | rel_close + + assert close.all() + + +def compare_in_memory( + process: pws.base.Process, + answers: dict[pws.base.adapter.AdapterNetcdf], + rtol: float = 1.0e-15, + atol: float = 1.0e-15, + equal_nan: bool = True, + strict: bool = False, + also_check_w_np: bool = True, + error_message: str = None, +): + # TODO: docstring + for var in process.get_variables(): + answers[var].advance() + assert_allclose( + process[var], + answers[var].current.data, + atol=atol, + rtol=rtol, + equal_nan=equal_nan, + strict=strict, + also_check_w_np=also_check_w_np, + error_message=error_message, + ) + + +def compare_netcdfs( + var_list: list, + results_dir: pl.Path, + answers_dir: pl.Path, + rtol: float = 1.0e-15, + atol: float = 1.0e-15, + equal_nan: bool = True, + strict: bool = False, + also_check_w_np: bool = True, + error_message: str = None, +): + # TODO: docstring + # TODO: improve error message + # TODO: collect failures in a try and report at end + for var in var_list: + answer = xr.open_dataarray(answers_dir / f"{var}.nc") + result = xr.open_dataarray(results_dir / f"{var}.nc") + + if error_message is None: + error_message = f"Comparison of variable '{var}' was unsuccessful" + + assert_allclose( + actual=result.values, + desired=answer.values, + rtol=rtol, + atol=atol, + equal_nan=equal_nan, + strict=strict, + also_check_w_np=also_check_w_np, + error_message=error_message, + ) diff --git a/bin/prms_linux_gfort_dbl_prec b/bin/prms_linux_gfort_dbl_prec new file mode 100755 index 00000000..00b2ab53 Binary files /dev/null and b/bin/prms_linux_gfort_dbl_prec differ diff --git a/bin/prms_linux b/bin/prms_linux_gfort_mixed_prec similarity index 100% rename from bin/prms_linux rename to bin/prms_linux_gfort_mixed_prec diff --git a/bin/prms_mac_intel_gfort_dbl_prec b/bin/prms_mac_intel_gfort_dbl_prec new file mode 100755 index 00000000..9abb32c4 Binary files /dev/null and b/bin/prms_mac_intel_gfort_dbl_prec differ diff --git a/bin/prms_mac b/bin/prms_mac_intel_gfort_mixed_prec similarity index 100% rename from bin/prms_mac rename to bin/prms_mac_intel_gfort_mixed_prec diff --git a/bin/prms_mac_m1_ifort_dbl_prec b/bin/prms_mac_m1_ifort_dbl_prec new file mode 100755 index 00000000..a9e1ee79 Binary files /dev/null and b/bin/prms_mac_m1_ifort_dbl_prec differ diff --git a/bin/prms_mac_m1_intel b/bin/prms_mac_m1_ifort_mixed_prec similarity index 99% rename from bin/prms_mac_m1_intel rename to bin/prms_mac_m1_ifort_mixed_prec index 7485b20d..56403fae 100755 Binary files a/bin/prms_mac_m1_intel and b/bin/prms_mac_m1_ifort_mixed_prec differ diff --git a/bin/prms_win.exe b/bin/prms_win.exe deleted file mode 100755 index dfd2476e..00000000 Binary files a/bin/prms_win.exe and /dev/null differ diff --git a/bin/prms_win_gfort_dbl_prec.exe b/bin/prms_win_gfort_dbl_prec.exe new file mode 100644 index 00000000..9ac58623 Binary files /dev/null and b/bin/prms_win_gfort_dbl_prec.exe differ diff --git a/examples/runoff_errors.ipynb b/examples/runoff_errors.ipynb new file mode 100644 index 00000000..c49d5172 --- /dev/null +++ b/examples/runoff_errors.ipynb @@ -0,0 +1,824 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4bf60102-e977-4b2c-82b1-86780c283c71", + "metadata": {}, + "source": [ + "# Runoff mass balance errors\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6936fb48-fb5c-463e-8642-11287a578ab1", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "# auto-format the code in this notebook\n", + "%load_ext jupyter_black" + ] + }, + { + "cell_type": "markdown", + "id": "c7b84519-933d-4284-838d-f5a8d4ec01f7", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d316a5de-19ef-4008-8c68-1e192c3f89d9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "import pathlib as pl\n", + "from pprint import pprint\n", + "from shutil import rmtree, copy2\n", + "\n", + "import hvplot.xarray # noqa\n", + "from icecream import ic\n", + "from IPython.display import display\n", + "import numpy as np\n", + "import pywatershed as pws\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c188b5b-fc91-480c-a1a2-ce1121f92915", + "metadata": {}, + "outputs": [], + "source": [ + "# Configuration\n", + "domain_name = \"drb_2yr\"\n", + "\n", + "nb_output_dir = pl.Path(\"./runoff_errors\")\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9e616a3-03a3-4a13-9a33-181af59f1c71", + "metadata": {}, + "outputs": [], + "source": [ + "pws_root = pws.constants.__pywatershed_root__\n", + "domain_dir = pws_root / f\"../test_data/{domain_name}\"\n", + "nb_output_dir.mkdir(exist_ok=True)\n", + "\n", + "zero = np.zeros([1])[0]\n", + "epsilon64 = np.finfo(zero).eps\n", + "epsilon32 = np.finfo(zero.astype(\"float32\")).eps" + ] + }, + { + "cell_type": "markdown", + "id": "0a05c33a-c86b-433a-a400-54c7f2df8eb9", + "metadata": {}, + "source": [ + "## Run PRMS double precision runs and convert to netcdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "daa94738-1c8b-4113-a529-c6dd0b4a0284", + "metadata": {}, + "outputs": [], + "source": [ + "bin_dir = pws_root / \"../bin/\"\n", + "# bin_mixed = bin_dir / \"prms_mac_m1_ifort_mixed_prec\"\n", + "bin_double = bin_dir / \"prms_mac_m1_ifort_dbl_prec\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "611c0d7a-7bfd-422d-ae53-3978b0c2eb83", + "metadata": {}, + "outputs": [], + "source": [ + "def run_prms(binary: pl.Path, run_dir: pl.Path, skip_if_exists=False):\n", + " import shlex\n", + " import subprocess\n", + "\n", + " from pywatershed import CsvFile, Soltab\n", + " from pywatershed.parameters import PrmsParameters\n", + "\n", + " if skip_if_exists and run_dir.exists():\n", + " print(\n", + " f\"Run ({run_dir}) already exists and skip_if_exists=True. Using existing run.\"\n", + " )\n", + " return None\n", + " \n", + " run_dir.mkdir() # must not exist, on user to delete\n", + " copy2(binary, run_dir / binary.name)\n", + " for ff in [\n", + " \"control.test\",\n", + " \"myparam.param\",\n", + " \"tmax.cbh\",\n", + " \"tmin.cbh\",\n", + " \"prcp.cbh\",\n", + " \"sf_data\",\n", + " ]:\n", + " copy2(domain_dir / ff, run_dir / ff)\n", + "\n", + " output_dir = run_dir / \"output\"\n", + " output_dir.mkdir()\n", + "\n", + " exe_command = f\"time ./{binary.name} control.test -MAXDATALNLEN 60000 2>&1 | tee run.log\"\n", + " result = subprocess.run(\n", + " exe_command,\n", + " shell=True,\n", + " # stdout = subprocess.PIPE,\n", + " stderr=subprocess.STDOUT,\n", + " universal_newlines=True,\n", + " cwd=run_dir,\n", + " )\n", + "\n", + " # these will be useful in what follows \n", + " params = pws.parameters.PrmsParameters.load(\n", + " domain_dir / \"myparam.param\"\n", + " ).parameters\n", + " \n", + " # convert to netcdf\n", + " # could make these arguments\n", + " chunking = {\n", + " \"time\": 0,\n", + " \"doy\": 0,\n", + " \"nhm_id\": 0,\n", + " \"nhm_seg\": 0,\n", + " }\n", + "\n", + " output_csvs = output_dir.glob(\"*.csv\")\n", + " for cc in output_csvs:\n", + " if cc.name in [\"stats.csv\"]:\n", + " continue\n", + " nc_path = cc.with_suffix(\".nc\")\n", + " CsvFile(cc).to_netcdf(nc_path, chunk_sizes=chunking)\n", + "\n", + " # previous and change variables\n", + " for vv in [\n", + " \"pk_ice\",\n", + " \"freeh2o\",\n", + " \"soil_moist\",\n", + " \"hru_impervstor\",\n", + " \"dprst_stor_hru\",\n", + " ]:\n", + " data = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n", + " prev_da = data.copy()\n", + " prev_da[:] = np.roll(prev_da.values, 1, axis=0)\n", + " assert (prev_da[1:, :].values == data[0:-1, :].values).all()\n", + " prev_da[0, :] = np.zeros(1)[\n", + " 0\n", + " ] # np.nan better but causes plotting to fail\n", + " change_da = data - prev_da\n", + " prev_da.rename(f\"{vv}_prev\").to_dataset().to_netcdf(\n", + " output_dir / f\"{vv}_prev.nc\"\n", + " )\n", + " data[f\"{vv}_prev\"] = xr.open_dataset(output_dir / f\"{vv}_prev.nc\")[\n", + " f\"{vv}_prev\"\n", + " ]\n", + "\n", + " change_da.rename(f\"{vv}_change\").to_dataset().to_netcdf(\n", + " output_dir / f\"{vv}_change.nc\"\n", + " )\n", + " data[f\"{vv}_change\"] = xr.open_dataset(output_dir / f\"{vv}_change.nc\")[\n", + " f\"{vv}_change\"\n", + " ]\n", + "\n", + " # through_rain\n", + " dep_vars = [\n", + " \"net_ppt\",\n", + " \"pptmix_nopack\",\n", + " \"snowmelt\",\n", + " \"pkwater_equiv\",\n", + " \"snow_evap\",\n", + " \"net_snow\",\n", + " \"net_rain\",\n", + " ]\n", + " data = {}\n", + " for vv in dep_vars:\n", + " data[vv] = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n", + "\n", + " nearzero = 1.0e-6\n", + "\n", + " cond1 = data[\"net_ppt\"] > zero\n", + " cond2 = data[\"pptmix_nopack\"] != 0\n", + " cond3 = data[\"snowmelt\"] < nearzero\n", + " cond4 = data[\"pkwater_equiv\"] < epsilon32\n", + " cond5 = data[\"snow_evap\"] < nearzero\n", + " cond6 = data[\"net_snow\"] < nearzero\n", + "\n", + " through_rain = data[\"net_rain\"] * zero\n", + " # these are in reverse order\n", + " through_rain[:] = np.where(\n", + " cond1 & cond3 & cond4 & cond6, data[\"net_rain\"], zero\n", + " )\n", + " through_rain[:] = np.where(\n", + " cond1 & cond3 & cond4 & cond5, data[\"net_ppt\"], through_rain\n", + " )\n", + " through_rain[:] = np.where(cond1 & cond2, data[\"net_rain\"], through_rain)\n", + "\n", + " through_rain.to_dataset(name=\"through_rain\").to_netcdf(\n", + " output_dir / \"through_rain.nc\"\n", + " )\n", + " through_rain.close()\n", + "\n", + " # infil_hru\n", + " imperv_frac = params[\"hru_percent_imperv\"]\n", + " dprst_frac = params[\"dprst_frac\"]\n", + " perv_frac = 1.0 - imperv_frac - dprst_frac\n", + " da = xr.open_dataset(output_dir / \"infil.nc\")[\"infil\"].rename(\"infil_hru\")\n", + " da *= perv_frac\n", + " da.to_dataset().to_netcdf(output_dir / \"infil_hru.nc\")\n", + " da.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e591db15-82cc-4843-a46b-87d037823700", + "metadata": {}, + "outputs": [], + "source": [ + "# run_prms(\n", + "# bin_mixed,\n", + "# nb_output_dir / f\"{domain_name}_prms_mixed_run\",\n", + "# skip_if_exists=skip_if_exists_prms_mixed,\n", + "# )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "791f0f08-06fb-445f-b6bb-9262b9f9a6f0", + "metadata": {}, + "outputs": [], + "source": [ + "# %debug" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "155fc8ac-da82-4df7-a33a-a744818d46d7", + "metadata": {}, + "outputs": [], + "source": [ + "prms_dbl_run_dir = nb_output_dir / f\"{domain_name}_prms_double_run\"\n", + "run_prms(\n", + " bin_double, prms_dbl_run_dir, skip_if_exists=skip_if_exists_prms_double\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "34bcdb09-c0c2-45c4-ab5f-fc22f94e77ef", + "metadata": {}, + "source": [ + "## Run pywatershed run forced with output from PRMS double precision run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01d0f43f-30cd-480e-9200-7c9507d4b263", + "metadata": {}, + "outputs": [], + "source": [ + "process = [pws.PRMSRunoff]\n", + "\n", + "pws_run_dir = nb_output_dir / f\"{domain_name}_pws_run\"\n", + "input_dir = pws_run_dir / \"pws_input\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80bc300c-2779-497f-877f-6a502531ca4a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "control = pws.Control.load(domain_dir / \"control.test\")\n", + "output_dir = pws_run_dir / \"output\"\n", + "control.options = control.options | {\n", + " \"input_dir\": input_dir,\n", + " \"budget_type\": \"error\",\n", + " \"calc_method\": \"numpy\",\n", + " \"netcdf_output_dir\": output_dir,\n", + "}\n", + "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a0312cf-a286-4247-b4f0-ca8f424a09d2", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "if output_dir.exists() and skip_if_exists_pws:\n", + " print(\n", + " f\"Output ({output_dir}) already exists and skip_if_exists=True. Using existing run.\"\n", + " )\n", + "\n", + "else:\n", + " input_dir.mkdir(exist_ok=True, parents=True)\n", + " for ff in prms_dbl_run_dir.glob(\"*.nc\"):\n", + " copy2(ff, input_dir / ff.name)\n", + " for ff in (prms_dbl_run_dir / \"output\").glob(\"*.nc\"):\n", + " copy2(ff, input_dir / ff.name)\n", + "\n", + " submodel = pws.Model(\n", + " process,\n", + " control=control,\n", + " parameters=params,\n", + " )\n", + "\n", + " submodel.run(finalize=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36d3958b-7c37-4e3e-9b46-cb29b9f75f0d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0e7d6ae-6502-4815-b534-1bcf9c2b948a", + "metadata": {}, + "outputs": [], + "source": [ + "for vv in process[0].get_variables():\n", + " print(vv)\n", + " assert (output_dir / f\"{vv}.nc\").exists()\n", + " try:\n", + " assert (input_dir / f\"{vv}.nc\").exists()\n", + " except:\n", + " print(f\"********** {vv} not in input_dir\")" + ] + }, + { + "cell_type": "markdown", + "id": "fa4ca04f-c639-479c-8832-8cea51fc79be", + "metadata": {}, + "source": [ + "## Start by comparing the budget variables" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f149cac-3343-4425-bddd-9e641245f3c9", + "metadata": {}, + "outputs": [], + "source": [ + "budget_terms = process[0].get_mass_budget_terms()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2cb23a3-38c3-488b-ace3-118d5fd1ae64", + "metadata": {}, + "outputs": [], + "source": [ + "# additional variables\n", + "budget_terms[\"outputs\"] += [\n", + " \"dprst_insroff_hru\",\n", + " \"dprst_stor_hru\",\n", + " \"contrib_fraction\",\n", + " \"infil\",\n", + " \"infil_hru\",\n", + " \"sroff\",\n", + " \"hru_sroffp\",\n", + " \"hru_sroffi\",\n", + " # \"imperv_stor\",\n", + " # \"imperv_evap\",\n", + " \"hru_impervevap\",\n", + " \"hru_impervstor\",\n", + " # \"hru_impervstor_old\",\n", + " \"hru_impervstor_change\",\n", + " # \"dprst_vol_frac\",\n", + " # \"dprst_vol_clos\",\n", + " # \"dprst_vol_open\",\n", + " # \"dprst_vol_clos_frac\",\n", + " # \"dprst_vol_open_frac\",\n", + " # \"dprst_area_clos\",\n", + " # \"dprst_area_open\",\n", + " # \"dprst_area_clos_max\",\n", + " # \"dprst_area_open_max\",\n", + " \"dprst_sroff_hru\",\n", + " \"dprst_seep_hru\",\n", + " \"dprst_evap_hru\",\n", + " \"dprst_insroff_hru\",\n", + " \"dprst_stor_hru\",\n", + " # \"dprst_stor_hru_old\",\n", + " \"dprst_stor_hru_change\",\n", + "]\n", + "budget_terms[\"outputs\"] = list(set(budget_terms[\"outputs\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b823fdf-e855-4263-b694-621c6b74005a", + "metadata": {}, + "outputs": [], + "source": [ + "comparisons = {}\n", + "for term, vars in budget_terms.items():\n", + " if term == \"inputs\":\n", + " continue\n", + " print(term)\n", + " for vv in vars:\n", + " print(\" \", vv)\n", + "\n", + " pws_file = output_dir / f\"{vv}.nc\"\n", + " assert (pws_file).exists()\n", + " pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n", + "\n", + " prms_file = input_dir / f\"{vv}.nc\"\n", + " assert (prms_file).exists()\n", + " prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n", + "\n", + " comparisons[vv] = xr.merge([pws_ds, prms_ds])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec5e0551-9d22-4522-b80e-4c2f54eedd28", + "metadata": {}, + "outputs": [], + "source": [ + "# comparisons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66f80562-8679-4ef3-b007-470692de064a", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_var(var_name, diff=False, nhm_id: list = None):\n", + " from textwrap import fill\n", + "\n", + " # lines = textwrap.wrap(text, width, break_long_words=False)\n", + " meta = pws.meta.find_variables(var_name)[var_name]\n", + " ylabel = f\"{fill(meta['desc'], 40)}\\n({meta['units']})\"\n", + " title = var_name\n", + " ds = comparisons[var_name]\n", + "\n", + " if diff:\n", + " ds = ds.copy()\n", + " ds[\"error\"] = ds[\"pws\"] - ds[\"prms\"]\n", + " ds[\"relative_error\"] = ds[\"error\"] / ds[\"prms\"]\n", + " del ds[\"pws\"], ds[\"prms\"]\n", + " ylabel = \"Difference PWS - PRMS\\n\" + ylabel\n", + " title = \"ERRORS: Difference in \" + title\n", + "\n", + " if (nhm_id is not None) and (len(nhm_id) > 0):\n", + " ds = ds.where(ds.nhm_id.isin(nhm_id), drop=True)\n", + "\n", + " display(\n", + " ds.hvplot(\n", + " frame_width=700,\n", + " groupby=\"nhm_id\",\n", + " title=title,\n", + " ylabel=ylabel,\n", + " # fontsize={\"ylabel\": \"9px\"},\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38a84f40-5684-4c73-a36f-027e3dc1bdbb", + "metadata": {}, + "outputs": [], + "source": [ + "def var_close(var_name):\n", + " print(var_name)\n", + " var_ds = comparisons[var_name]\n", + " abs_diff = abs(var_ds[\"pws\"] - var_ds[\"prms\"])\n", + " rel_abs_diff = abs_diff / var_ds[\"prms\"]\n", + " rtol = atol = 1.0e-7\n", + " close = (abs_diff < atol) | (rel_abs_diff < rtol)\n", + " if close.all():\n", + " plot_var(var_name, diff=False)\n", + "\n", + " else:\n", + " wh_not_close = np.where(~close)\n", + " nhm_ids = abs_diff.nhm_id[wh_not_close[1]]\n", + " plot_var(var_name, diff=True, nhm_id=nhm_ids)\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68c77573-6c17-45a1-8730-6cba7e832958", + "metadata": {}, + "outputs": [], + "source": [ + "var_close(\"hru_impervstor_change\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3b9d5ff-05e1-4468-bef6-54b3da7f5b0c", + "metadata": {}, + "outputs": [], + "source": [ + "for var_name in comparisons.keys():\n", + " var_close(var_name)" + ] + }, + { + "cell_type": "markdown", + "id": "62a7b9e8-5e95-4aa6-867c-1ad48f04b5a6", + "metadata": {}, + "source": [ + "## Look at specific time and budget errors\n", + "\n", + "> UserWarning: The flux unit balance not equal to the change in unit storageat time 1979-04-26T00:00:00 \n", + "> and at the following locations for PRMSRunoff: (array([341, 509]),)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ce47dd4-fc1d-4ba3-ba59-b69d8c692409", + "metadata": {}, + "outputs": [], + "source": [ + "budget_terms = process[0].get_mass_budget_terms()\n", + "budget_terms[\"inputs\"] += [\n", + " \"net_ppt\",\n", + " \"net_rain\",\n", + " \"net_snow\",\n", + " \"pptmix_nopack\",\n", + " \"pk_ice_prev\",\n", + " \"freeh2o_prev\",\n", + " \"newsnow\",\n", + " \"snow_evap\",\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f52f1fb-585e-49bb-9db9-b2e87d62a973", + "metadata": {}, + "outputs": [], + "source": [ + "budget_cases = [\n", + " (\n", + " \"1979-01-11T00:00:00\",\n", + " [0, 1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13],\n", + " ),\n", + "]\n", + "\n", + "case_ind = 0\n", + "budget_time = np.datetime64(budget_cases[case_ind][0])\n", + "budget_location_inds = budget_cases[case_ind][1]\n", + "\n", + "budget_comps = {}\n", + "for term, vars in budget_terms.items():\n", + " print(term)\n", + "\n", + " for vv in vars:\n", + " print(\" \", vv)\n", + "\n", + " if term == \"inputs\":\n", + " pws_file = input_dir / f\"{vv}.nc\"\n", + " else:\n", + " pws_file = output_dir / f\"{vv}.nc\"\n", + "\n", + " assert (pws_file).exists()\n", + " pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n", + "\n", + " prms_file = input_dir / f\"{vv}.nc\"\n", + " assert (prms_file).exists()\n", + " prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n", + "\n", + " budget_comps[vv] = (\n", + " xr.merge([pws_ds, prms_ds])\n", + " .sel(time=budget_time)\n", + " .isel(nhm_id=budget_location_inds)\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0dad7836-6157-4982-a77e-f0f46210a27f", + "metadata": {}, + "outputs": [], + "source": [ + "bc = budget_comps\n", + "inputs = bc[\"through_rain\"] + bc[\"snowmelt\"] + bc[\"intcp_changeover\"]\n", + "outputs = (\n", + " bc[\"hru_sroffi\"]\n", + " + bc[\"hru_sroffp\"]\n", + " + bc[\"dprst_sroff_hru\"]\n", + " + bc[\"infil_hru\"]\n", + " + bc[\"hru_impervevap\"]\n", + " + bc[\"dprst_seep_hru\"]\n", + " + bc[\"dprst_evap_hru\"]\n", + ")\n", + "storage_changes = bc[\"hru_impervstor_change\"] + bc[\"dprst_stor_hru_change\"]\n", + "balance = inputs - outputs - storage_changes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "127606d6-f28f-43c1-90f0-4f3c9bad3a5f", + "metadata": {}, + "outputs": [], + "source": [ + "ic(budget_location_inds)\n", + "ic(inputs.prms.values)\n", + "ic(outputs.prms.values)\n", + "ic(storage_changes.prms.values)\n", + "\n", + "ic(\"-----------\")\n", + "\n", + "ic(bc[\"through_rain\"].pws.values)\n", + "\n", + "ic(bc[\"snow_evap\"].prms.values)\n", + "ic(bc[\"hru_impervstor_change\"].prms.values)\n", + "ic(bc[\"hru_impervstor_change\"].pws.values)\n", + "ic(bc[\"dprst_stor_hru_change\"].prms.values)\n", + "ic(bc[\"dprst_stor_hru_change\"].pws.values)\n", + "ic(balance.prms.values)\n", + "\n", + "# ic(bc[\"hru_sroffi\"].prms.sum().values)\n", + "# ic(bc[\"hru_sroffp\"].prms.sum().values)\n", + "# ic(bc[\"dprst_sroff_hru\"].prms.sum().values)\n", + "# ic(bc[\"infil_hru\"].prms.sum().values)\n", + "# ic(bc[\"hru_impervevap\"].prms.sum().values)\n", + "# ic(bc[\"dprst_seep_hru\"].prms.sum().values)\n", + "# ic(bc[\"dprst_evap_hru\"].prms.sum().values)\n", + "\n", + "# ic(storage_changes.prms.values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63297394-ec68-4645-babb-a46cd6b680bf", + "metadata": {}, + "outputs": [], + "source": [ + "ic((balance - bc[\"through_rain\"]).pws.values)\n", + "ic((balance - bc[\"through_rain\"]).prms.values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52d6a2ed-d900-4f66-bf0b-e75c1d11fce9", + "metadata": {}, + "outputs": [], + "source": [ + "((balance.pws - balance.prms) < 1.0e-8).all().values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3b27604-a410-4e8c-ac77-a136aac554df", + "metadata": {}, + "outputs": [], + "source": [ + "balance.pws.sum()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f993bb1-7590-4e47-a2a7-b0bc10350ef3", + "metadata": {}, + "outputs": [], + "source": [ + "ic(bc[\"through_rain\"].pws.values)\n", + "ic(bc[\"net_rain\"].pws.values)\n", + "ic(bc[\"net_snow\"].pws.values)\n", + "ic(bc[\"net_ppt\"].pws.values)\n", + "ic(bc[\"pptmix_nopack\"].pws.values)\n", + "ic(bc[\"newsnow\"].pws.values)\n", + "ic((bc[\"pk_ice_prev\"].pws.values + bc[\"freeh2o_prev\"].pws.values) < epsilon32)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "968ffe3e-a984-4b36-bc9e-7f6a360b245c", + "metadata": {}, + "outputs": [], + "source": [ + "input_max = max(\n", + " abs(bc[\"through_rain\"]), abs(bc[\"snowmelt\"]), +abs(bc[\"intcp_changeover\"])\n", + ")\n", + "output_max = max(\n", + " abs(bc[\"hru_sroffi\"]),\n", + " abs(bc[\"hru_sroffp\"]),\n", + " abs(bc[\"dprst_sroff_hru\"]),\n", + " abs(bc[\"infil_hru\"]),\n", + " abs(bc[\"hru_impervevap\"]),\n", + " abs(bc[\"dprst_seep_hru\"]),\n", + " abs(bc[\"dprst_evap_hru\"]),\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "203b6241-601f-4a4d-b2d7-d0358435be9c", + "metadata": {}, + "outputs": [], + "source": [ + "input_max.pws.values.tolist()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bfbe75c-a67e-4d7a-b29d-3c48d04cc53b", + "metadata": {}, + "outputs": [], + "source": [ + "output_max.prms.values.tolist()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a0e4a51-6794-4aae-b379-c658f7c73bce", + "metadata": {}, + "outputs": [], + "source": [ + "print((balance.pws / output_max.pws.max()).values)\n", + "print((balance.prms / output_max.prms.max()).values)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1e738e0-2761-4d07-8b86-a2c8a5403c5c", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython" + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/snow_errors.ipynb b/examples/snow_errors.ipynb new file mode 100644 index 00000000..b8243a21 --- /dev/null +++ b/examples/snow_errors.ipynb @@ -0,0 +1,773 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4bf60102-e977-4b2c-82b1-86780c283c71", + "metadata": {}, + "source": [ + "# Snow mass balance errors\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6936fb48-fb5c-463e-8642-11287a578ab1", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "# auto-format the code in this notebook\n", + "%load_ext jupyter_black" + ] + }, + { + "cell_type": "markdown", + "id": "c7b84519-933d-4284-838d-f5a8d4ec01f7", + "metadata": {}, + "source": [ + "## Setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d316a5de-19ef-4008-8c68-1e192c3f89d9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "import pathlib as pl\n", + "from pprint import pprint\n", + "from shutil import rmtree, copy2\n", + "\n", + "import hvplot.xarray # noqa\n", + "from icecream import ic\n", + "from IPython.display import display\n", + "import numpy as np\n", + "import pywatershed as pws\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c188b5b-fc91-480c-a1a2-ce1121f92915", + "metadata": {}, + "outputs": [], + "source": [ + "domain_name = \"ucb_2yr\"\n", + "\n", + "nb_output_dir = pl.Path(\"./snow_errors\")\n", + "\n", + "skip_if_exists_prms_mixed = True\n", + "skip_if_exists_prms_double = True\n", + "skip_if_exists_pws = True" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f807cf6-f982-4b5a-91eb-89294bf1c08e", + "metadata": {}, + "outputs": [], + "source": [ + "pws_root = pws.constants.__pywatershed_root__\n", + "domain_dir = pws_root / f\"../test_data/{domain_name}\"\n", + "nb_output_dir.mkdir(exist_ok=True)\n", + "\n", + "zero = pws.constants.zero\n", + "epsilon64 = pws.constants.epsilon64\n", + "epsilon32 = pws.constants.epsilon32" + ] + }, + { + "cell_type": "markdown", + "id": "0a05c33a-c86b-433a-a400-54c7f2df8eb9", + "metadata": {}, + "source": [ + "## Run PRMS mixed and double precision runs and convert to netcdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "daa94738-1c8b-4113-a529-c6dd0b4a0284", + "metadata": {}, + "outputs": [], + "source": [ + "bin_dir = pws_root / \"../bin/\"\n", + "# bin_mixed = bin_dir / \"prms_521_mixed_mac_m1_intel\"\n", + "bin_double = bin_dir / \"prms_mac_m1_ifort_dbl_prec\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "611c0d7a-7bfd-422d-ae53-3978b0c2eb83", + "metadata": {}, + "outputs": [], + "source": [ + "def run_prms(binary: pl.Path, run_dir: pl.Path, skip_if_exists=False):\n", + " import shlex\n", + " import subprocess\n", + "\n", + " from pywatershed import CsvFile, Soltab\n", + "\n", + " from pywatershed.parameters import PrmsParameters\n", + "\n", + " if skip_if_exists and run_dir.exists():\n", + " print(\n", + " f\"Run ({run_dir}) already exists and skip_if_exists=True. Using existing run.\"\n", + " )\n", + " return None\n", + "\n", + " run_dir.mkdir() # must not exist, on user to delete\n", + " copy2(binary, run_dir / binary.name)\n", + " for ff in [\n", + " \"control.test\",\n", + " \"myparam.param\",\n", + " \"tmax.cbh\",\n", + " \"tmin.cbh\",\n", + " \"prcp.cbh\",\n", + " \"sf_data\",\n", + " ]:\n", + " copy2(domain_dir / ff, run_dir / ff)\n", + "\n", + " output_dir = run_dir / \"output\"\n", + " output_dir.mkdir()\n", + "\n", + " exe_command = f\"time ./{binary.name} control.test -MAXDATALNLEN 60000 2>&1 | tee run.log\"\n", + " result = subprocess.run(\n", + " exe_command,\n", + " shell=True,\n", + " # stdout = subprocess.PIPE,\n", + " stderr=subprocess.STDOUT,\n", + " universal_newlines=True,\n", + " cwd=run_dir,\n", + " )\n", + "\n", + " # these will be useful in what follows\n", + " params = pws.parameters.PrmsParameters.load(\n", + " domain_dir / \"myparam.param\"\n", + " ).parameters\n", + "\n", + " # convert to netcdf\n", + " # could make these arguments\n", + " chunking = {\n", + " \"time\": 0,\n", + " \"doy\": 0,\n", + " \"nhm_id\": 100,\n", + " \"nhm_seg\": 100,\n", + " }\n", + "\n", + " output_csvs = output_dir.glob(\"*.csv\")\n", + " for cc in output_csvs:\n", + " if cc.name in [\"stats.csv\"]:\n", + " continue\n", + " nc_path = cc.with_suffix(\".nc\")\n", + " CsvFile(cc).to_netcdf(nc_path, chunk_sizes=chunking)\n", + "\n", + " # solar tables\n", + " soltab_file = run_dir / \"soltab_debug\"\n", + " # the nhm_ids are not available in the solta_debug file currently, so get\n", + " # them from the domain parameters\n", + " params = PrmsParameters.load(run_dir / \"myparam.param\")\n", + " nhm_ids = params.parameters[\"nhm_id\"]\n", + "\n", + " soltab = Soltab(\n", + " soltab_file,\n", + " output_dir=output_dir,\n", + " nhm_ids=nhm_ids,\n", + " chunk_sizes=chunking,\n", + " )\n", + "\n", + " for var in soltab.variables:\n", + " assert (output_dir / f\"{var}.nc\").exists()\n", + "\n", + " # previous and change variables\n", + " for vv in [\"pk_ice\", \"freeh2o\", \"soil_moist\"]:\n", + " data = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n", + " prev_da = data.copy()\n", + " prev_da[:] = np.roll(prev_da.values, 1, axis=0)\n", + " assert (prev_da[1:, :].values == data[0:-1, :].values).all()\n", + " prev_da[0, :] = np.zeros(1)[\n", + " 0\n", + " ] # np.nan better but causes plotting to fail\n", + " change_da = data - prev_da\n", + " prev_da.rename(f\"{vv}_prev\").to_dataset().to_netcdf(\n", + " output_dir / f\"{vv}_prev.nc\"\n", + " )\n", + " data[f\"{vv}_prev\"] = xr.open_dataset(output_dir / f\"{vv}_prev.nc\")[\n", + " f\"{vv}_prev\"\n", + " ]\n", + "\n", + " change_da.rename(f\"{vv}_change\").to_dataset().to_netcdf(\n", + " output_dir / f\"{vv}_change.nc\"\n", + " )\n", + " data[f\"{vv}_change\"] = xr.open_dataset(output_dir / f\"{vv}_change.nc\")[\n", + " f\"{vv}_change\"\n", + " ]\n", + " # through_rain\n", + " dep_vars = [\n", + " \"net_ppt\",\n", + " \"pptmix_nopack\",\n", + " \"snowmelt\",\n", + " \"pkwater_equiv\",\n", + " \"snow_evap\",\n", + " \"net_snow\",\n", + " \"net_rain\",\n", + " ]\n", + " data = {}\n", + " for vv in dep_vars:\n", + " data[vv] = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n", + "\n", + " nearzero = 1.0e-6\n", + "\n", + " cond1 = data[\"net_ppt\"] > zero\n", + " cond2 = data[\"pptmix_nopack\"] != 0\n", + " cond3 = data[\"snowmelt\"] < nearzero\n", + " cond4 = data[\"pkwater_equiv\"] < epsilon32\n", + " cond5 = data[\"snow_evap\"] < nearzero\n", + " cond6 = data[\"net_snow\"] < nearzero\n", + "\n", + " through_rain = data[\"net_rain\"] * zero\n", + " # these are in reverse order\n", + " through_rain[:] = np.where(\n", + " cond1 & cond3 & cond4 & cond6, data[\"net_rain\"], zero\n", + " )\n", + " through_rain[:] = np.where(\n", + " cond1 & cond3 & cond4 & cond5, data[\"net_ppt\"], through_rain\n", + " )\n", + " through_rain[:] = np.where(cond1 & cond2, data[\"net_rain\"], through_rain)\n", + "\n", + " through_rain.to_dataset(name=\"through_rain\").to_netcdf(\n", + " output_dir / \"through_rain.nc\"\n", + " )\n", + " through_rain.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e591db15-82cc-4843-a46b-87d037823700", + "metadata": {}, + "outputs": [], + "source": [ + "# run_prms(bin_mixed, nb_output_dir / \"prms_mixed_run\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "155fc8ac-da82-4df7-a33a-a744818d46d7", + "metadata": {}, + "outputs": [], + "source": [ + "prms_dbl_run_dir = nb_output_dir / f\"{domain_name}_prms_double_run\"\n", + "run_prms(\n", + " bin_double, prms_dbl_run_dir, skip_if_exists=skip_if_exists_prms_double\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "34bcdb09-c0c2-45c4-ab5f-fc22f94e77ef", + "metadata": {}, + "source": [ + "## pywatershed run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01d0f43f-30cd-480e-9200-7c9507d4b263", + "metadata": {}, + "outputs": [], + "source": [ + "process = [pws.PRMSSnow]\n", + "pws_run_dir = nb_output_dir / f\"{domain_name}_pws_run\"\n", + "input_dir = pws_run_dir / \"pws_input\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80bc300c-2779-497f-877f-6a502531ca4a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + } + }, + "outputs": [], + "source": [ + "control = pws.Control.load(domain_dir / \"control.test\")\n", + "output_dir = pws_run_dir / \"output\"\n", + "control.options = control.options | {\n", + " \"input_dir\": input_dir,\n", + " \"budget_type\": \"warn\",\n", + " \"calc_method\": \"numpy\",\n", + " \"netcdf_output_dir\": output_dir,\n", + "}\n", + "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65276857-f940-46df-8385-a9791bc79894", + "metadata": {}, + "outputs": [], + "source": [ + "if output_dir.exists() and skip_if_exists_pws:\n", + " print(\n", + " f\"Output ({output_dir}) already exists and skip_if_exists=True. Using existing run.\"\n", + " )\n", + "\n", + "else:\n", + " input_dir.mkdir(exist_ok=True, parents=True)\n", + " for ff in prms_dbl_run_dir.glob(\"*.nc\"):\n", + " copy2(ff, input_dir / ff.name)\n", + " for ff in (prms_dbl_run_dir / \"output\").glob(\"*.nc\"):\n", + " copy2(ff, input_dir / ff.name)\n", + "\n", + " submodel = pws.Model(\n", + " process,\n", + " control=control,\n", + " parameters=params,\n", + " )\n", + "\n", + " submodel.run(finalize=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a056b6a-d31a-492b-a761-3a33f42ba08c", + "metadata": {}, + "outputs": [], + "source": [ + "%debug" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0e7d6ae-6502-4815-b534-1bcf9c2b948a", + "metadata": {}, + "outputs": [], + "source": [ + "for vv in process[0].get_variables():\n", + " print(vv)\n", + " assert (output_dir / f\"{vv}.nc\").exists()\n", + " try:\n", + " assert (input_dir / f\"{vv}.nc\").exists()\n", + " except:\n", + " print(f\"********** {vv} not in input_dir\")" + ] + }, + { + "cell_type": "markdown", + "id": "fa4ca04f-c639-479c-8832-8cea51fc79be", + "metadata": {}, + "source": [ + "## Start by comparing the budget variables" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f149cac-3343-4425-bddd-9e641245f3c9", + "metadata": {}, + "outputs": [], + "source": [ + "budget_terms = pws.PRMSSnow.get_mass_budget_terms()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2cb23a3-38c3-488b-ace3-118d5fd1ae64", + "metadata": {}, + "outputs": [], + "source": [ + "# additional variables\n", + "budget_terms[\"outputs\"] += [\n", + " \"pk_ice_prev\",\n", + " \"freeh2o_prev\",\n", + " \"newsnow\",\n", + " \"pptmix_nopack\",\n", + " # \"ai\",\n", + " \"albedo\",\n", + " # 'frac_swe',\n", + " \"freeh2o\",\n", + " \"freeh2o_change\",\n", + " \"freeh2o_prev\",\n", + " #' iasw',\n", + " # 'int_alb',\n", + " \"iso\",\n", + " # 'lso',\n", + " # 'lst',\n", + " # \"mso\",\n", + " \"newsnow\",\n", + " \"pk_def\",\n", + " \"pk_den\",\n", + " \"pk_depth\",\n", + " \"pk_ice\",\n", + " \"pk_ice_change\",\n", + " \"pk_ice_prev\",\n", + " # 'pk_precip',\n", + " \"pk_temp\",\n", + " # 'pksv',\n", + " \"pkwater_ante\",\n", + " \"pkwater_equiv\",\n", + " # 'pkwater_equiv_change',\n", + " \"pptmix_nopack\",\n", + " # 'pss',\n", + " \"pst\",\n", + " # \"salb\",\n", + " #' scrv',\n", + " #' slst',\n", + " \"snow_evap\",\n", + " \"snowcov_area\",\n", + " # 'snowcov_areasv',\n", + " \"snowmelt\",\n", + " # 'snsv',\n", + " \"tcal\",\n", + " \"through_rain\",\n", + "]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b823fdf-e855-4263-b694-621c6b74005a", + "metadata": {}, + "outputs": [], + "source": [ + "comparisons = {}\n", + "for term, vars in budget_terms.items():\n", + " if term == \"inputs\":\n", + " continue\n", + " print(term)\n", + " for vv in sorted(vars):\n", + " print(\" \", vv)\n", + "\n", + " pws_file = output_dir / f\"{vv}.nc\"\n", + " assert (pws_file).exists()\n", + " pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n", + "\n", + " prms_file = input_dir / f\"{vv}.nc\"\n", + " assert prms_file.exists()\n", + " prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n", + "\n", + " comparisons[vv] = xr.merge([pws_ds, prms_ds])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec5e0551-9d22-4522-b80e-4c2f54eedd28", + "metadata": {}, + "outputs": [], + "source": [ + "# comparisons" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66f80562-8679-4ef3-b007-470692de064a", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_var(var_name, diff=False, nhm_id: list = None):\n", + " from textwrap import fill\n", + "\n", + " # lines = textwrap.wrap(text, width, break_long_words=False)\n", + " meta = pws.meta.find_variables(var_name)[var_name]\n", + " ylabel = f\"{fill(meta['desc'], 40)}\\n({meta['units']})\"\n", + " title = var_name\n", + " ds = comparisons[var_name]\n", + "\n", + " if diff:\n", + " ds = ds.copy()\n", + " ds[\"error\"] = ds[\"pws\"] - ds[\"prms\"]\n", + " # ds[\"relative_error\"] = ds[\"error\"] / ds[\"prms\"]\n", + " # ds[\"relative_error\"] = xr.where(\n", + " # abs(ds[\"prms\"]) < 1.0e-7, np.nan, ds[\"relative_error\"]\n", + " # )\n", + " del ds[\"pws\"], ds[\"prms\"]\n", + " ylabel = \"Difference PWS - PRMS\\n\" + ylabel\n", + " title = \"ERRORS: Difference in \" + title\n", + "\n", + " if (nhm_id is not None) and (len(nhm_id) > 0):\n", + " ds = ds.where(ds.nhm_id.isin(nhm_id), drop=True)\n", + "\n", + " display(\n", + " ds.hvplot(\n", + " frame_width=700,\n", + " groupby=\"nhm_id\",\n", + " title=title,\n", + " ylabel=ylabel,\n", + " # fontsize={\"ylabel\": \"9px\"},\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38a84f40-5684-4c73-a36f-027e3dc1bdbb", + "metadata": {}, + "outputs": [], + "source": [ + "def var_close(var_name):\n", + " print(var_name)\n", + " var_ds = comparisons[var_name]\n", + " abs_diff = abs(var_ds[\"pws\"] - var_ds[\"prms\"])\n", + " rel_abs_diff = abs_diff / var_ds[\"prms\"]\n", + " rtol = 0.02\n", + " atol = 1.0e-2\n", + " close = (abs_diff < atol) | (rel_abs_diff < rtol)\n", + " if close.all():\n", + " plot_var(var_name, diff=False)\n", + "\n", + " else:\n", + " wh_not_close = np.where(~close)\n", + " nhm_ids = abs_diff.nhm_id[wh_not_close[1]]\n", + " plot_var(var_name, diff=False, nhm_id=nhm_ids)\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3b9d5ff-05e1-4468-bef6-54b3da7f5b0c", + "metadata": {}, + "outputs": [], + "source": [ + "if False:\n", + " for var_name in comparisons.keys():\n", + " var_close(var_name)" + ] + }, + { + "cell_type": "markdown", + "id": "da70848a-3187-4ae7-b045-125c83e451c9", + "metadata": {}, + "source": [ + "## Mass-balance errors\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5d5e323-a739-43a6-a372-5bcb7dc7b775", + "metadata": {}, + "outputs": [], + "source": [ + "budget_terms = pws.PRMSSnow.get_mass_budget_terms()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ff49bb08-0ae4-41e9-9e52-adeb0d48cbbe", + "metadata": {}, + "outputs": [], + "source": [ + "budget_terms" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bacee4a6-0952-4c96-aa50-b05f5664e74b", + "metadata": {}, + "outputs": [], + "source": [ + "prms_mass = {}\n", + "for term, vars in budget_terms.items():\n", + " prms_mass[term] = {}\n", + " for vv in vars:\n", + " prms_mass[term][vv] = xr.open_dataarray(input_dir / f\"{vv}.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3deaaa12-3094-4283-8edf-1c1ad0acdfe7", + "metadata": {}, + "outputs": [], + "source": [ + "prms_mass_sums = {}\n", + "for term, vars in prms_mass.items():\n", + " prms_mass_sums[term] = None\n", + " for vv in vars:\n", + " data = xr.open_dataarray(input_dir / f\"{vv}.nc\")\n", + " if prms_mass_sums[term] is None:\n", + " prms_mass_sums[term] = data\n", + " else:\n", + " prms_mass_sums[term] += xr.open_dataarray(input_dir / f\"{vv}.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "552ff261-c959-4d89-8f7f-53160c928f75", + "metadata": {}, + "outputs": [], + "source": [ + "prms_mass_balance = (\n", + " prms_mass_sums[\"inputs\"]\n", + " - prms_mass_sums[\"outputs\"]\n", + " - prms_mass_sums[\"storage_changes\"]\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf5228fd-0a6c-4c9e-bc29-adf8373743f2", + "metadata": {}, + "outputs": [], + "source": [ + "prms_mass_balance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbd6c480-fc75-4546-9246-2c5f5b2aee46", + "metadata": {}, + "outputs": [], + "source": [ + "wh_imbalance = np.where(prms_mass_balance.values > 1e-5)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2be69429-2ab7-4307-90c6-dedb4262635a", + "metadata": {}, + "outputs": [], + "source": [ + "wh_imbalance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bb8a452-ae77-48e9-9120-02a8c1ebcd75", + "metadata": {}, + "outputs": [], + "source": [ + "prms_mass_balance.values[wh_imbalance]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c979a40-5740-4971-a0aa-d5e0f5dd88b7", + "metadata": {}, + "outputs": [], + "source": [ + "freeh2o = xr.open_dataarray(input_dir / \"freeh2o.nc\")\n", + "pk_ice = xr.open_dataarray(input_dir / \"pk_ice.nc\")\n", + "pkwater_equiv = xr.open_dataarray(input_dir / \"pkwater_equiv.nc\")\n", + "freeh2o_prev = xr.open_dataarray(input_dir / \"freeh2o_prev.nc\")\n", + "pk_ice_prev = xr.open_dataarray(input_dir / \"pk_ice_prev.nc\")\n", + "pkwater_ante = xr.open_dataarray(input_dir / \"pkwater_ante.nc\")\n", + "freeh2o_change = xr.open_dataarray(input_dir / \"freeh2o_change.nc\")\n", + "pk_ice_change = xr.open_dataarray(input_dir / \"pk_ice_change.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be45ba5c-75ce-49ac-89d1-8f91aaaadec3", + "metadata": {}, + "outputs": [], + "source": [ + "ic(pk_ice.values[wh_imbalance])\n", + "ic(pk_ice_prev.values[wh_imbalance])\n", + "ic(pk_ice_change.values[wh_imbalance])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3220c630-52c8-4074-8794-54a9f462a7a3", + "metadata": {}, + "outputs": [], + "source": [ + "ic(freeh2o.values[wh_imbalance])\n", + "ic(freeh2o_prev.values[wh_imbalance])\n", + "ic(freeh2o_change.values[wh_imbalance])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7ae9821-c6b6-422b-89d8-04f1037c83eb", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "np.set_printoptions(edgeitems=30, linewidth=100000)\n", + "\n", + "for term, vars in prms_mass.items():\n", + " print(term)\n", + " for vv in vars:\n", + " print(f\" {vv}\")\n", + " print(f\" {prms_mass[term][vv].values[wh_imbalance]}\")\n", + "print(\"imbalance\")\n", + "print(f\" {prms_mass_balance.values[wh_imbalance]}\")\n", + "print(\"freeh2o prev\")\n", + "print(f\" {freeh2o_prev.values[wh_imbalance]}\")\n", + "print(\"freeh2o\")\n", + "print(f\" {freeh2o.values[wh_imbalance]}\")\n", + "print(\"pk_ice prev\")\n", + "print(f\" {pk_ice_prev.values[wh_imbalance]}\")\n", + "print(\"pk_ice\")\n", + "print(f\" {pk_ice.values[wh_imbalance]}\")\n", + "print(\"pkwater_ante\")\n", + "print(f\" {pkwater_ante.values[wh_imbalance]}\")\n", + "print(\"pkwater_equiv\")\n", + "print(f\" {pkwater_equiv.values[wh_imbalance]}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a882b844-6fa3-4005-8660-36d45deb51f6", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython" + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/prms_src/prms5.2.1/makelist b/prms_src/prms5.2.1/makelist index 6f9ffa48..b4b0e1e6 100644 --- a/prms_src/prms5.2.1/makelist +++ b/prms_src/prms5.2.1/makelist @@ -37,7 +37,7 @@ endif $(info ----------------) $(info FC is $(FC)) -$(info FC is $(CC)) +$(info CC is $(CC)) $(info DBL_PREC is false) $(info ----------------) $(info ) diff --git a/prms_src/prms5.2.1/makelist_double_precision b/prms_src/prms5.2.1/makelist_double_precision index b45d4f4f..223f23cc 100644 --- a/prms_src/prms5.2.1/makelist_double_precision +++ b/prms_src/prms5.2.1/makelist_double_precision @@ -37,7 +37,7 @@ endif $(info ----------------) $(info FC is $(FC)) -$(info FC is $(CC)) +$(info CC is $(CC)) $(info DBL_PREC is true) $(info ----------------) $(info ) @@ -78,7 +78,7 @@ endif # This is a shortcut to alleviate errors when parsing parameter inputs in soltab. ifeq ($(detected_OS), Windows) ifeq ($(FC), gfortran) - FFLAGS = $(OPTLEVEL) -fno-second-underscore -fallow-argument-mismatch + FFLAGS = -freal-4-real-8 $(OPTLEVEL) -fno-second-underscore -fallow-argument-mismatch endif else ifeq ($(FC), gfortran) @@ -99,7 +99,7 @@ endif ifeq ($(detected_OS), Windows) ifeq ($(CC), gcc) - CFLAGS = $(OPTLEVEL) -D$(ARC) -D_UF -Wall + CFLAGS = -Dfloat=double $(OPTLEVEL) -D$(ARC) -D_UF -Wall endif else ifeq ($(CC), gcc) diff --git a/prms_src/prms5.2.1/prms/snowcomp.f90 b/prms_src/prms5.2.1/prms/snowcomp.f90 index fad3cf19..ef2bd151 100644 --- a/prms_src/prms5.2.1/prms/snowcomp.f90 +++ b/prms_src/prms5.2.1/prms/snowcomp.f90 @@ -1435,7 +1435,8 @@ INTEGER FUNCTION snorun() ENDIF ENDIF - ENDIF + ENDIF ! ends iff statement between steps 3 and 4 about 265 lines above + ! LAST check to clear out all arrays if packwater is gone IF ( Pkwater_equiv(i)<=0.0D0 ) THEN @@ -2598,6 +2599,9 @@ SUBROUTINE snowevap(Potet_sublim, Potet, Snowcov_area, Snow_evap, & !RAPCOMMENT - CHANGED TO CHECK FOR NEGATIVE PACK ICE ! If all pack ice is removed, then there cannot be a ! heat deficit + ! JLM mass balance fix, if we take pkwater_equiv - ez below + ! JLM we need to take ez from both pk_ice and freeh2o here + freeh2o = freeh2o + pk_ice Pk_ice = 0.0 Pk_def = 0.0 Pk_temp = 0.0 @@ -2613,17 +2617,22 @@ SUBROUTINE snowevap(Potet_sublim, Potet, Snowcov_area, Snow_evap, & Pkwater_equiv = Pkwater_equiv - ez Snow_evap = ez ENDIF + IF ( Snow_evap<0.0 ) THEN Pkwater_equiv = Pkwater_equiv - DBLE(Snow_evap) + IF ( Pkwater_equiv<0.0D0 ) THEN IF ( Print_debug>DEBUG_less ) THEN IF ( Pkwater_equiv<-DNEARZERO ) & & PRINT *, 'snowpack issue, negative pkwater_equiv in snowevap', Pkwater_equiv - Pkwater_equiv = 0.0D0 + Pkwater_equiv = 0.0D0 ! JLM: this is INSIDE a debug statement? will change the answers + ! is this in the originial source ENDIF ENDIF + Snow_evap = 0.0 ENDIF + avail_et = Potet - Hru_intcpevap - Snow_evap IF ( avail_et<0.0 ) THEN ! PRINT *, 'snow evap', snow_evap, avail_et, pkwater_equiv diff --git a/pywatershed/base/budget.py b/pywatershed/base/budget.py index a36e99d7..d8a428e2 100644 --- a/pywatershed/base/budget.py +++ b/pywatershed/base/budget.py @@ -8,7 +8,7 @@ from pywatershed.base.control import Control -from ..constants import zero +from ..constants import epsilon64, one, zero from ..utils.formatting import pretty_print from ..utils.netcdf_utils import NetCdfWrite from .accessor import Accessor @@ -298,18 +298,27 @@ def _sum_storage_changes(self): return self._sum("storage_changes") def _calc_unit_balance(self): - unit_balance = self._inputs_sum - self._outputs_sum self._zero_sum = True - if not np.allclose( - unit_balance, - self._storage_changes_sum, - rtol=self.rtol, - atol=self.atol, - ): + + # roll our own np.allclose so we can diagnose the not close points + unit_balance = self._inputs_sum - self._outputs_sum + abs_diff = np.abs(unit_balance - self._storage_changes_sum) + mask_div_zero = self._storage_changes_sum < epsilon64 + rel_abs_diff = abs_diff / np.where( + mask_div_zero, one, self._storage_changes_sum + ) + abs_close = abs_diff < self.atol + rel_close = rel_abs_diff < self.rtol + either_close = abs_close | rel_close + cond = np.where(mask_div_zero, abs_close, either_close) + + if not cond.all(): self._zero_sum = False + wh_not_cond = np.where(~cond) msg = ( - "The flux unit balance not equal to the change in unit storage" - f": {self.description}" + "The flux unit balance not equal to the change in unit " + f"storage at time {self.control.current_time} and at the " + f"following locations for {self.description}: {wh_not_cond}" ) if self.imbalance_fatal: raise ValueError(msg) diff --git a/pywatershed/base/control.py b/pywatershed/base/control.py index d99f80b4..3f3730cf 100644 --- a/pywatershed/base/control.py +++ b/pywatershed/base/control.py @@ -108,19 +108,19 @@ def current_year(self): return datetime_year(self._current_time) @property - def current_month(self): - """Get the current month.""" - return datetime_month(self._current_time) + def current_month(self, zero_based: bool = False): + """Get the current month in 1-12 (unless zero based).""" + return datetime_month(self._current_time, zero_based=zero_based) @property - def current_doy(self): - """Get the current day of year.""" - return datetime_doy(self._current_time) + def current_doy(self, zero_based: bool = False): + """Get the current day of year in 1-366 (unless zero based).""" + return datetime_doy(self._current_time, zero_based=zero_based) @property - def current_dowy(self): - """Get the current day of water year.""" - return datetime_dowy(self._current_time) + def current_dowy(self, zero_based: bool = False): + """Get the current day of water year in 1-366 (unless zero-based).""" + return datetime_dowy(self._current_time, zero_based=zero_based) @property def current_epiweek(self): diff --git a/pywatershed/hydrology/prms_runoff.py b/pywatershed/hydrology/prms_runoff.py index a26f9f06..85304316 100644 --- a/pywatershed/hydrology/prms_runoff.py +++ b/pywatershed/hydrology/prms_runoff.py @@ -7,7 +7,14 @@ from ..base.adapter import adaptable from ..base.conservative_process import ConservativeProcess from ..base.control import Control -from ..constants import HruType, nan, numba_num_threads, zero +from ..constants import ( + HruType, + epsilon32, + epsilon64, + nan, + numba_num_threads, + zero, +) from ..parameters import Parameters NEARZERO = 1.0e-6 @@ -25,6 +32,8 @@ LAND = HruType.LAND.value LAKE = HruType.LAKE.value +# TODO: using through_rain and not net_rain and net_ppt is a WIP + class PRMSRunoff(ConservativeProcess): """PRMS surface runoff. @@ -188,7 +197,7 @@ def get_init_values() -> dict: "hru_impervstor_change": zero, "dprst_vol_frac": zero, "dprst_vol_clos": zero, - "dprst_vol_open": zero, + "dprst_vol_open": nan, "dprst_vol_clos_frac": zero, "dprst_vol_open_frac": zero, "dprst_area_clos": zero, @@ -209,6 +218,7 @@ def get_mass_budget_terms(): return { "inputs": [ # "net_rain", + # "net_snow", "through_rain", "snowmelt", "intcp_changeover", @@ -500,6 +510,7 @@ def _calculate(self, time_length, vectorized=False): compute_infil=self.compute_infil, dprst_comp=self.dprst_comp, imperv_et=self.imperv_et, + through_rain=self.through_rain, ) self.infil_hru[:] = self.infil * self.hru_frac_perv @@ -586,10 +597,8 @@ def _calculate_numpy( compute_infil, dprst_comp, imperv_et, + through_rain, ): - # move towards replacing net_rain and other variables. - # self.net_rain[:] = self.through_rain - dprst_chk = 0 infil[:] = 0.0 for k in prange(nhru): @@ -615,8 +624,7 @@ def _calculate_numpy( hru_impervevap[i] = 0.0 avail_et = potet[i] - snow_evap[i] - hru_intcpevap[i] - # apparently not used - # availh2o = intcp_changeover[i] + net_rain[i] + availh2o = intcp_changeover[i] + net_rain[i] ( sri, @@ -648,6 +656,7 @@ def _calculate_numpy( srp=srp, check_capacity=check_capacity, perv_comp=perv_comp, + through_rain=through_rain[i], ) dprst_flag = ACTIVE # cdl todo: hardwired @@ -709,7 +718,7 @@ def _calculate_numpy( dprst_flow_coef=dprst_flow_coef[i], dprst_seep_rate_clos=dprst_seep_rate_clos[i], avail_et=avail_et, - net_rain=net_rain[i], + net_rain=availh2o, dprst_in=dprst_in[i], srp=srp, sri=sri, @@ -816,12 +825,14 @@ def compute_infil( srp, check_capacity, perv_comp, + through_rain, ): isglacier = False # todo -- hardwired - cascade_active = False hru_flag = 0 if hru_type == LAND or isglacier: hru_flag = 1 + + cascade_active = False if cascade_active: raise Exception("bad bad bad") else: @@ -846,17 +857,21 @@ def compute_infil( # if rain/snow event with no antecedent snowpack, # compute the runoff from the rain first and then proceed with the # snowmelt computations - if pptmix_nopack == ACTIVE: - avail_water = avail_water + net_rain - infil = infil + net_rain + # double_counting = 0 + cond2 = pptmix_nopack != 0 + # if pptmix_nopack == ACTIVE: + if cond2: + # double_counting += 1 + avail_water = avail_water + through_rain + infil = infil + through_rain if hru_flag == 1: infil, srp, contrib_fraction = perv_comp( soil_moist_prev, carea_max, smidx_coef, smidx_exp, - net_rain, - net_rain, + through_rain, + through_rain, infil, srp, ) @@ -866,13 +881,19 @@ def compute_infil( # procedure is used. If there is no snowpack and no precip, # then check for melt from last of snowpack. If rain/snow mix # with no antecedent snowpack, compute snowmelt portion of runoff. + cond3 = snowmelt < NEARZERO + cond4 = pkwater_equiv < epsilon64 + cond1 = net_ppt > zero + cond6 = net_snow < NEARZERO + if snowmelt > 0.0: + # if not cond3: # this results in less precision accuracy # There was no snowmelt but a snowpack may exist. If there is # no snowpack then check for rain on a snowfree HRU. avail_water = avail_water + snowmelt infil = infil + snowmelt if hru_flag == 1: - if pkwater_equiv > 0.0 or net_ppt - net_snow < NEARZERO: + if (pkwater_equiv > 0.0) or (net_rain < NEARZERO): # Pervious area computations infil, srp = check_capacity( soil_moist_prev, @@ -883,6 +904,11 @@ def compute_infil( ) else: # Snowmelt occurred and depleted the snowpack + # this frequently gets counted along with pptmix_nopack + # double_counting += 1 + # if double_counting > 1: + # print("snowmelt") + infil, srp, contrib_fraction = perv_comp( soil_moist_prev, carea_max, @@ -894,20 +920,27 @@ def compute_infil( srp, ) - elif pkwater_equiv < DNEARZERO: + # elif pkwater_equiv < DNEARZERO: + elif cond4: # If no snowmelt and no snowpack but there was net snow then # snowpack was small and was lost to sublimation. - if net_snow < NEARZERO and net_rain > 0.0: - avail_water = avail_water + net_rain - infil = infil + net_rain + # if net_snow < NEARZERO and net_rain > 0.0: + if cond6 and through_rain > 0.0: + # cond3 & cond4 & cond6 & cond1 + # this is through_rain's top/most narrow case + avail_water = avail_water + through_rain + infil = infil + through_rain + # double_counting += 1 + # if double_counting > 1: + # print("cond4") if hru_flag == 1: infil, srp, contrib_fraction = perv_comp( soil_moist_prev, carea_max, smidx_coef, smidx_exp, - net_rain, - net_rain, + through_rain, + through_rain, infil, srp, ) @@ -925,6 +958,9 @@ def compute_infil( srp, ) + # if double_counting > 1: + # assert False + if hruarea_imperv > 0.0: imperv_stor = imperv_stor + avail_water if hru_flag == 1: @@ -1048,7 +1084,9 @@ def dprst_comp( if sri < 0.0: if sri < -NEARZERO: sri = 0.0 - dprst_insroff_hru = dprst_srp + dprst_sri + + # <<< + dprst_insroff_hru = dprst_srp + dprst_sri dprst_area_open = 0.0 if dprst_vol_open > 0.0: diff --git a/pywatershed/hydrology/prms_snow.py b/pywatershed/hydrology/prms_snow.py index 86694cde..f46f970d 100644 --- a/pywatershed/hydrology/prms_snow.py +++ b/pywatershed/hydrology/prms_snow.py @@ -64,6 +64,7 @@ maxalb = 15 ONETHIRD = 1.0 / 3.0 +tiny_snowpack = 1.0e-12 tcind = 0 @@ -697,22 +698,24 @@ def _calculate_numpy( # cals = zero # JLM this is unnecessary. - # newsnow - newsnow[:] = False - net_snow_gt_zero = net_snow > zero - wh_net_snow_gt_zero = np.where(net_snow_gt_zero) - newsnow[wh_net_snow_gt_zero] = True + # newsnow is a doganostic for prms_snow, so it lives here + newsnow = np.where(net_snow > zero, True, False) - # default assumption - pptmix_nopack[:] = False + # JLM TODO: there's a conditional here we dont have + # in fotran trd is scalar and the RHS terms are vector? + trd = orad_hru / soltab_horad_potsw + + frac_swe[:] = zero + pk_precip[:] = zero # [inches] + snowmelt[:] = zero # [inches] + snow_evap[:] = zero # [inches] + tcal[:] = zero + ai[:] = zero for jj in prange(nhru): if hru_type[jj] == HruType.LAKE.value: continue - # < - trd = orad_hru[jj] / soltab_horad_potsw[jj] - # If it's the first julian day of the water year, several # variables need to be reset: # - reset the previous snow water eqivalent plus new snow to 0 @@ -733,15 +736,8 @@ def _calculate_numpy( # TIME PERIOD # ************************************************************** - # By default, the precipitation added to snowpack, snowmelt, - # and snow evaporation are 0. - # JLM: this could happen outside the loop - pk_precip[jj] = zero # [inches] - snowmelt[jj] = zero # [inches] - snow_evap[jj] = zero # [inches] - frac_swe[jj] = zero - ai[jj] = zero - tcal[jj] = zero + # default assumption + pptmix_nopack[jj] = False # If the day of the water year is beyond the forced melt day # indicated by the parameter, then set the flag indicating melt @@ -764,20 +760,15 @@ def _calculate_numpy( # print(f"pk_ice 0 : {pk_ice[dbgind]}") # print(f"tcal 0 : {tcal[dbgind]}") - if pkwater_equiv[jj] < epsilon64: - # No existing snowpack - if not newsnow[jj]: - # Skip the HRU if there is no snowpack and no new snow - # Reset to be sure it is zero if snowpack melted on last - # timestep. - snowcov_area[jj] = zero - continue - else: - # We ahave new snow; the initial snow-covered area is - # complete (1) - # JLM: why set this here? just for the case of no existing - # snow? This might be removable. - snowcov_area[jj] = one # [fraction of area] + if (pkwater_equiv[jj] < epsilon64) and (newsnow[jj] == 0): + # Skip the HRU if there is no snowpack and no new snow + # Reset to be sure it is zero if snowpack melted on last + # timestep. + snowcov_area[jj] = zero + continue + + if newsnow[jj] and (pkwater_equiv[jj] < epsilon64): + snowcov_area[jj] = one # << # HRU STEP 1 - DEAL WITH PRECIPITATION AND ITS EFFECT ON THE WATER @@ -785,6 +776,7 @@ def _calculate_numpy( # **************************************************************** month_ind = current_month - 1 + # PRMS conditonal moved inside function ( freeh2o[jj], iasw[jj], @@ -906,13 +898,16 @@ def _calculate_numpy( # ) # print(f"tcal 3 : {tcal[dbgind]}") + if pkwater_equiv[jj] > zero: # HRU STEP 4 - DETERMINE RADIATION FLUXES AND SNOWPACK # STATES NECESSARY FOR ENERGY BALANCE # ********************************************************** ( freeh2o[jj], + iasw[jj], iso[jj], lso[jj], + mso[jj], pk_def[jj], pk_den[jj], pk_depth[jj], @@ -921,8 +916,9 @@ def _calculate_numpy( pkwater_equiv[jj], pss[jj], tcal[jj], + snowmelt[jj], ) = calc_step_4( - trd, + trd[jj], calc_calin=calc_calin, calc_caloss=calc_caloss, calc_snowbal=calc_snowbal, @@ -1035,13 +1031,13 @@ def _calculate_numpy( # insufficient to reset albedo, then reduce the cumulative # new snow by the amount melted during the period (but # don't let it be negative). - if lst[jj]: + if lst[jj] > 0: snsv[jj] = snsv[jj] - snowmelt[jj] if snsv[jj] < zero: snsv[jj] = zero - # <<<< + # <<<< The if starting between steps 3 & 4 ends here # LAST check to clear out all arrays if packwater is gone if pkwater_equiv[jj] <= zero: # if verbose: @@ -1073,17 +1069,30 @@ def _calculate_numpy( snowcov_areasv[jj] = zero # rsr, not in original code ai[jj] = zero frac_swe[jj] = zero + scrv[jj] = zero + pksv[jj] = zero - # << + # << end of space loop and previous if pkwater_equiv_change[:] = pkwater_equiv - pkwater_ante freeh2o_change[:] = freeh2o - freeh2o_prev pk_ice_change[:] = pk_ice - pk_ice_prev - wh_through = ( - ((pk_ice_prev + freeh2o_prev) <= epsilon64) & ~newsnow - ) | (pptmix_nopack == 1) - through_rain[:] = np.where(wh_through, net_rain, zero) + nearzero = 1.0e-6 + cond1 = net_ppt > zero + cond2 = pptmix_nopack != 0 + cond3 = snowmelt < nearzero + cond4 = pkwater_equiv < epsilon32 + cond5 = snow_evap < nearzero + cond6 = net_snow < nearzero + # reverse order from the if statements + through_rain[:] = np.where( + cond1 & cond3 & cond4 & cond6, net_rain, zero + ) + through_rain[:] = np.where( + cond1 & cond3 & cond4 & cond5, net_ppt, through_rain + ) + through_rain[:] = np.where(cond1 & cond2, net_rain, through_rain) return ( ai, @@ -1203,10 +1212,10 @@ def _calc_ppt_to_pack( snowmelt, ) - caln = zero - calpr = zero - calps = zero - pndz = zero + # caln = zero + # calpr = zero + # calps = zero + # pndz = zero tsnow = tavgc # [degrees C] @@ -1575,16 +1584,6 @@ def _calc_calin( pk_def = pk_def - cal # [cal/cm^2] pk_temp = -1 * pk_def / (pkwater_equiv * 1.27) # [degrees C] - elif abs(dif) < epsilon32: - # JLM: I moved this from an else at the bottom to this - # point in the conditional chain. - # JLM: The test had been equality with zero, changed to - # less than epsilon32. - # (2) Just enough heat to overcome heat deficit - # Set temperature and heat deficit to zero. the pack is "ripe" - pk_temp = zero # [degrees C] - pk_def = zero # [cal/cm^2] - elif dif > zero: # (3) More than enough heat to overcome heat deficit (melt ice)... # Calculate the potential amount of snowmelt from excess heat in @@ -1701,6 +1700,17 @@ def _calc_calin( pss = pkwater_equiv # [inches] # <<< + else: # abs(dif) < epsilon64: + # JLM: The test had been equality with zero, changed to + # less than epsilon64. + # (2) Just enough heat to overcome heat deficit + # Set temperature and heat deficit to zero. the pack is "ripe" + pk_temp = zero # [degrees C] + pk_def = zero # [cal/cm^2] + + if not (pkwater_equiv > zero): + pk_den = zero + return ( freeh2o, iasw, @@ -1854,10 +1864,11 @@ def _calc_snowcov( # < # Calculate the ratio of the current packwater equivalent to the # maximum packwater equivalent for the given snowpack. - if ai == zero: - frac_swe = zero - else: + if ai > epsilon64: frac_swe = pkwater_equiv / ai # [fraction] + frac_swe = min(one, frac_swe) + else: + frac_swe = zero # < # There are 3 potential conditions for the snow area curve: @@ -2204,7 +2215,6 @@ def _calc_snalbedo( # reset the albedo. # Reset the albedo states. slst = zero # [days] - # lst = 0 # [flag] lst = False # [flag] snsv = zero # [inches] @@ -2235,7 +2245,7 @@ def _calc_snalbedo( # threshold. # 4 options below (if-then, elseif, elseif, else) - if pptmix < one: + if pptmix == 0: # (3.1) This is not a mixed event... # During the accumulation season, the threshold for resetting # the albedo does not apply if there is a snow-only event. @@ -2439,12 +2449,6 @@ def _calc_step_4( # Save the current value of emissivity esv = emis # [fraction of radiation] - # The incoming shortwave radiation is the HRU radiation adjusted by the - # albedo (some is reflected back into the atmoshphere) and the - # transmission coefficient (some is intercepted by the winter - # vegetative canopy) - swn = swrad * (one - albedo) * rad_trncf # [cal/cm^2] or [Langleys] - # Set the convection-condensation for a half-day interval cec = cecn_coef * 0.5 # [cal/(cm^2 degC)] or [Langleys/degC] @@ -2454,49 +2458,6 @@ def _calc_step_4( # RSR: cov_type==4 is valid for trees (coniferous) cec = cec * 0.5 # [cal/(cm^2 degC)] or [Langleys/degC] - # < - # Calculate the new snow depth (Riley et al. 1973) - # RSR: the following 3 lines of code were developed by Rob Payn, - # on 7/10/2013 - # The snow depth depends on the previous snow pack water equivalent - # plus the current net snow. - pss = pss + net_snow # [inches] - dpt_before_settle = pk_depth + net_snow * deninv - dpt1 = dpt_before_settle + settle_const * ( - (pss * denmaxinv) - dpt_before_settle - ) - - # RAPCOMMENT - CHANGED TO THE APPROPRIATE FINITE DIFFERENCE - # APPROXIMATION OF SNOW DEPTH - # JLM: pk_depth is prognostic here - pk_depth = dpt1 # [inches] - - # Calculate the snowpack density - if dpt1 > zero: - pk_den = pkwater_equiv / dpt1 - else: - pk_den = zero # [inch water equiv / inch depth] - - # < - # The effective thermal conductivity is approximated (empirically) - # as zero077 times (snowpack density)^2 [cal / (sec g degC)] Therefore, - # the effective conductivity term (inside the square root) in the - # equation for conductive heat exchange can be calculated as follows: - # (zero077 * pk_den^2) / (pk_den * 0.5) - # where 0.5 is the specific heat of ice [cal / (g degC)] - # this simplifies to the following - effk = 0.0154 * pk_den # [unitless] - - # 13751 is the number of seconds in 12 hours over pi - # So for a half day, to calculate the conductive heat exchange per cm - # of snow per cm^2 area per degree temperature difference is the - # following - # In effect, multiplying cst times the temperature gradient gives the - # heatexchange by heat conducted (calories) per square cm of snowpack - cst = pk_den * ( - np.sqrt(effk * 13751.0) - ) # [cal/(cm^2 degC)] or [Langleys / degC] - # Check whether to force spring melt # Spring melt is forced if time is before the melt-force day and after # the melt-look day (parameters). If between these dates, the spring @@ -2539,75 +2500,139 @@ def _calc_step_4( # Set the flag indicating night time niteda = 1 # [flag] - # No shortwave (solar) radiation at night. - sw = zero # [cal / cm^2] or [Langleys] - # Temperature is halfway between the minimum and average temperature # for the day. temp = (tminc + tavgc) * 0.5 - # Track total heat flux from both night and day periods + if pkwater_equiv > zero: + # The incoming shortwave radiation is the HRU radiation adjusted by + # the albedo (some is reflected back into the atmoshphere) and the + # transmission coefficient (some is intercepted by the winter + # vegetative canopy) + swn = ( + swrad * (one - albedo) * rad_trncf + ) # [cal/cm^2] or [Langleys] + + # Calculate the new snow depth (Riley et al. 1973) + # RSR: the following 3 lines of code were developed by Rob Payn, + # on 7/10/2013 + # The snow depth depends on the previous snow pack water equivalent + # plus the current net snow. + pss = pss + net_snow # [inches] + # deninv is the inverse of den_init + dpt_before_settle = pk_depth + net_snow * deninv + dpt1 = dpt_before_settle + settle_const * ( + (pss * denmaxinv) - dpt_before_settle + ) - # Calculate the night time energy balance - ( - tcal, - freeh2o, - pk_def, - pk_ice, - pk_temp, - pkwater_equiv, - ) = calc_snowbal( - niteda=niteda, - cec=cec, - cst=cst, - esv=esv, - sw=sw, - temp=temp, - trd=trd, - calc_calin=calc_calin, - calc_caloss=calc_caloss, - canopy_covden=canopy_covden, - den_max=den_max, - denmaxinv=denmaxinv, - emis_noppt=emis_noppt, - freeh2o=freeh2o, - freeh2o_cap=freeh2o_cap, - hru_ppt=hru_ppt, - iasw=iasw, - pk_def=pk_def, - pk_den=pk_den, - pk_depth=pk_depth, - pk_ice=pk_ice, - pk_temp=pk_temp, - pkwater_equiv=pkwater_equiv, - pss=pss, - pst=pst, - snowcov_area=snowcov_area, - snowmelt=snowmelt, - tcal=tcal, - tstorm_mo=tstorm_mo, - ) + # RAPCOMMENT - CHANGED TO THE APPROPRIATE FINITE DIFFERENCE + # APPROXIMATION OF SNOW DEPTH + # JLM: pk_depth is prognostic here + pk_depth = dpt1 # [inches] + + # Calculate the snowpack density + if dpt1 > zero: + pk_den = pkwater_equiv / dpt1 + else: + pk_den = zero # [inch water equiv / inch depth] + + # < + # The effective thermal conductivity is approximated (empirically) + # as zero077 times (snowpack density)^2 [cal / (sec g degC)] Therefore, + # the effective conductivity term (inside the square root) in the + # equation for conductive heat exchange can be calculated as follows: + # (zero077 * pk_den^2) / (pk_den * 0.5) + # where 0.5 is the specific heat of ice [cal / (g degC)] + # this simplifies to the following + effk = 0.0154 * pk_den # [unitless] + + # 13751 is the number of seconds in 12 hours over pi + # So for a half day, to calculate the conductive heat exchange per cm + # of snow per cm^2 area per degree temperature difference is the + # following + # In effect, multiplying cst times the temperature gradient gives the + # heatexchange by heat conducted (calories) per square cm of snowpack + cst = pk_den * ( + np.sqrt(effk * 13751.0) + ) # [cal/(cm^2 degC)] or [Langleys / degC] + + # No shortwave (solar) radiation at night. + sw = zero # [cal / cm^2] or [Langleys] + + # Calculate the night time energy balance + ( + tcal, + freeh2o, + iasw, + pk_def, + pk_den, + pk_ice, + pk_depth, + pk_temp, + pss, + pst, + snowmelt, + pkwater_equiv, + ) = calc_snowbal( + niteda=niteda, + cec=cec, + cst=cst, + esv=esv, + sw=sw, + temp=temp, + trd=trd, + calc_calin=calc_calin, + calc_caloss=calc_caloss, + canopy_covden=canopy_covden, + den_max=den_max, + denmaxinv=denmaxinv, + emis_noppt=emis_noppt, + freeh2o=freeh2o, + freeh2o_cap=freeh2o_cap, + hru_ppt=hru_ppt, + iasw=iasw, + pk_def=pk_def, + pk_den=pk_den, + pk_depth=pk_depth, + pk_ice=pk_ice, + pk_temp=pk_temp, + pkwater_equiv=pkwater_equiv, + pss=pss, + pst=pst, + snowcov_area=snowcov_area, + snowmelt=snowmelt, + tcal=tcal, + tstorm_mo=tstorm_mo, + ) + # tcal set directly by calc_snowbal above [cal/cm^2] or [Langleys] - # [cal/cm^2] or [Langleys] + # < + # iswn = zero on ly a glacier variable # Compute energy balance for day period (if the snowpack still exists) # THIS SHOULD HAPPEN IN SNOBAL - if pkwater_equiv > zero: - # Set the flag indicating daytime - niteda = 2 # [flag] + # Set the flag indicating daytime + niteda = 2 # [flag] + # Temperature is halfway between the maximum and average + # temperature for the day. + temp = (tmaxc + tavgc) * 0.5 # [degrees C] + if pkwater_equiv > zero: # Set shortwave radiation as calculated earlier sw = swn # [cal/cm^2] or [Langleys] - # Temperature is halfway between the maximum and average - # temperature for the day. - temp = (tmaxc + tavgc) * 0.5 # [degrees C] ( cals, freeh2o, + iasw, pk_def, + pk_den, pk_ice, + pk_depth, pk_temp, + pss, + pst, + snowmelt, pkwater_equiv, ) = calc_snowbal( niteda=niteda, @@ -2647,8 +2672,10 @@ def _calc_step_4( # < return ( freeh2o, + iasw, iso, lso, + mso, pk_def, pk_den, pk_depth, @@ -2657,6 +2684,7 @@ def _calc_step_4( pkwater_equiv, pss, tcal, + snowmelt, ) @staticmethod @@ -2832,9 +2860,15 @@ def _calc_snowbal( return ( cal, freeh2o, + iasw, pk_def, + pk_den, pk_ice, + pk_depth, pk_temp, + pss, + pst, + snowmelt, pkwater_equiv, ) @@ -3003,9 +3037,15 @@ def _calc_snowbal( return ( cal, freeh2o, + iasw, pk_def, + pk_den, pk_ice, + pk_depth, pk_temp, + pss, + pst, + snowmelt, pkwater_equiv, ) @@ -3048,8 +3088,6 @@ def _calc_snowevap( # Set the evaporation to the pack water equivalent and set all # snowpack # variables to no-snowpack values. - snow_evap = pkwater_equiv # [inches] - snow_evap = pkwater_equiv # [inches] pkwater_equiv = zero # [inches] pk_ice = zero # [inches] @@ -3068,6 +3106,10 @@ def _calc_snowevap( # RAPCOMMENT - CHANGED TO CHECK FOR NEGATIVE PACK ICE # If all pack ice is removed, then there cannot be a heat # deficit. + + # JLM: mass balance fix only in our 5.2.1 prms version + freeh2o = freeh2o + pk_ice + pk_ice = zero pk_def = zero pk_temp = zero @@ -3095,8 +3137,10 @@ def _calc_snowevap( # f"snowevap: {pkwater_equiv}" # ) # # << - - pkwater_equiv = zero + # zero of pkwater_equiv is inside a debug statement that + # is dosent seem to be triggered in test runs + # pkwater_equiv = zero + pass # < snow_evap = zero @@ -3114,12 +3158,11 @@ def _calc_snowevap( # if verbose: # if pkwater_equiv < -epsilon64: # print( - # "snowpack issue 2, negative pkwater_equiv in " + # "snowpack issue 2, negative pkwater_equiv in" # f"snowevap: {pkwater_equiv}" # ) # # << - - # To be sure negative snowpack is ignored + # This zero of pkwater IS NOT in the debug statement pkwater_equiv = zero # < diff --git a/pywatershed/utils/csv_utils.py b/pywatershed/utils/csv_utils.py index 73a8dcf7..a0c9fffc 100644 --- a/pywatershed/utils/csv_utils.py +++ b/pywatershed/utils/csv_utils.py @@ -124,7 +124,7 @@ def to_netcdf( clobber: bool = True, zlib: bool = True, complevel: int = 4, - chunk_sizes: dict = {"time": 30, "hruid": 0}, + chunk_sizes: dict = {"time": 30, "nhm_id": 0, "nhm_seg": 0}, ) -> None: """Output the csv output data to a netcdf file @@ -192,14 +192,17 @@ def to_netcdf( ids = self._coordinates[dim_name] nids = len(ids) + dims = ("time", dim_name) + chunk_sizes_var = [chunk_sizes[vv] for vv in dims] + var = ds.createVariable( variable_name, variable_type, - ("time", dim_name), + dims, fill_value=nc4.default_fillvals[variable_type], # JLM: sus zlib=zlib, complevel=complevel, - chunksizes=tuple(chunk_sizes.values()), + chunksizes=tuple(chunk_sizes_var), ) # add additional meta data if self.meta.is_available(variable_name): diff --git a/pywatershed/utils/prms5util.py b/pywatershed/utils/prms5util.py index b3133125..a732651b 100644 --- a/pywatershed/utils/prms5util.py +++ b/pywatershed/utils/prms5util.py @@ -249,6 +249,7 @@ def __init__( soltab_file: fileish, output_dir: fileish = None, nhm_ids: np.ndarray = None, + chunk_sizes: dict = {"doy": 0, "nhm_id": 0}, ): self.soltab_file = soltab_file self.output_dir = output_dir @@ -273,7 +274,7 @@ def __init__( ) = load_soltab_debug(self.soltab_file) if self.output_dir: - self.to_netcdf() + self.to_netcdf(chunk_sizes=chunk_sizes) return @@ -283,7 +284,7 @@ def to_netcdf( nhm_ids: np.ndarray = None, zlib: bool = True, complevel: int = 4, - chunk_sizes: dict = {"doy": 0, "hruid": 0}, + chunk_sizes: dict = {"doy": 0, "nhm_id": 0}, ): # This is just different enough to make it it's own thing compared # to the CSV outputs of PRMS. @@ -320,14 +321,17 @@ def to_netcdf( ) hruid[:] = self.spatial_coord + dims = ("doy", self.spatial_coord_name) + chunk_sizes_var = [chunk_sizes[vv] for vv in dims] + variables[vv] = ds.createVariable( vv, "f4", - ("doy", self.spatial_coord_name), + dims, fill_value=nc4.default_fillvals["f4"], zlib=zlib, complevel=complevel, - chunksizes=tuple(chunk_sizes.values()), + chunksizes=tuple(chunk_sizes_var), ) variables[vv][:] = self[vv] ds.close() diff --git a/pywatershed/utils/time_utils.py b/pywatershed/utils/time_utils.py index 87581b90..dda65c45 100644 --- a/pywatershed/utils/time_utils.py +++ b/pywatershed/utils/time_utils.py @@ -18,27 +18,34 @@ def datetime_year(dt64: np.datetime64) -> int: return dt64.astype("datetime64[Y]").astype(int) + 1970 -def datetime_month(dt64: np.datetime64) -> int: +def datetime_month(dt64: np.datetime64, zero_based=False) -> int: """Get the month from np.datetime64""" - return dt64.astype("datetime64[M]").astype(int) % 12 + 1 + return dt64.astype("datetime64[M]").astype(int) % 12 + _offset(zero_based) -def datetime_doy(dt64: np.datetime64) -> int: +def datetime_doy(dt64: np.datetime64, zero_based=False) -> int: """Get day of year from np.datetime64""" return (dt64 - dt64.astype("datetime64[Y]")).astype( "timedelta64[D]" - ).astype(int) + 1 + ).astype(int) + _offset(zero_based) -def datetime_dowy(dt64: np.datetime64) -> int: +def datetime_dowy(dt64: np.datetime64, zero_based=False) -> int: """Get day of water year from np.datetime64""" year_start = datetime_year(dt64) if datetime_month(dt64) < 10: year_start -= 1 diff = dt64 - np.datetime64(f"{year_start}-10-01") - return diff.astype("timedelta64[D]").astype(int) + return diff.astype("timedelta64[D]").astype(int) + _offset(zero_based) def datetime_epiweek(dt64: np.datetime64) -> int: """Get CDC eipweek [1, 53] from np.datetime64""" return ew.Week.fromdate(dt64_to_dt(dt64)).week + + +def _offset(zero_based: bool): + if zero_based: + return 0 + else: + return 1 diff --git a/test_data/common/control.multi_hru b/test_data/common/control.multi_hru index df624841..effd9667 100644 --- a/test_data/common/control.multi_hru +++ b/test_data/common/control.multi_hru @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -78 +83 #### nhruOut_freq 1 @@ -251,13 +251,14 @@ nhruOut_format 4 #### nhruOutVar_names -78 +83 4 albedo cap_infil_tot contrib_fraction dprst_seep_hru dprst_evap_hru +dprst_insroff_hru dprst_sroff_hru dprst_stor_hru freeh2o @@ -288,7 +289,11 @@ net_snow newsnow orad_hru perv_actet +pk_def +pk_den +pk_depth pk_ice +pk_temp pkwater_ante pkwater_equiv potet @@ -345,7 +350,7 @@ nsegmentOutON_OFF nsegmentOutVars 1 1 -9 +10 #### nsegmentOut_freq 1 @@ -353,8 +358,9 @@ nsegmentOut_freq 1 #### nsegmentOutVar_names -9 +10 4 +seg_inflow seg_outflow seg_gwflow seg_sroff diff --git a/test_data/common/control.single_hru b/test_data/common/control.single_hru index e692a8b6..74ea87b3 100644 --- a/test_data/common/control.single_hru +++ b/test_data/common/control.single_hru @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -76 +81 #### nhruOut_freq 1 @@ -251,12 +251,13 @@ nhruOut_format 5 #### nhruOutVar_names -76 +81 4 albedo cap_infil_tot contrib_fraction dprst_evap_hru +dprst_insroff_hru dprst_seep_hru dprst_stor_hru dprst_sroff_hru @@ -288,7 +289,11 @@ net_snow newsnow orad_hru perv_actet +pk_def +pk_den +pk_depth pk_ice +pk_temp pkwater_ante pkwater_equiv potet @@ -343,7 +348,7 @@ nsegmentOutON_OFF nsegmentOutVars 1 1 -9 +10 #### nsegmentOut_freq 1 @@ -351,8 +356,9 @@ nsegmentOut_freq 1 #### nsegmentOutVar_names -9 +10 4 +seg_inflow seg_outflow seg_gwflow seg_sroff diff --git a/test_data/drb_2yr/control.test b/test_data/drb_2yr/control.test index df624841..effd9667 100644 --- a/test_data/drb_2yr/control.test +++ b/test_data/drb_2yr/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -78 +83 #### nhruOut_freq 1 @@ -251,13 +251,14 @@ nhruOut_format 4 #### nhruOutVar_names -78 +83 4 albedo cap_infil_tot contrib_fraction dprst_seep_hru dprst_evap_hru +dprst_insroff_hru dprst_sroff_hru dprst_stor_hru freeh2o @@ -288,7 +289,11 @@ net_snow newsnow orad_hru perv_actet +pk_def +pk_den +pk_depth pk_ice +pk_temp pkwater_ante pkwater_equiv potet @@ -345,7 +350,7 @@ nsegmentOutON_OFF nsegmentOutVars 1 1 -9 +10 #### nsegmentOut_freq 1 @@ -353,8 +358,9 @@ nsegmentOut_freq 1 #### nsegmentOutVar_names -9 +10 4 +seg_inflow seg_outflow seg_gwflow seg_sroff diff --git a/test_data/scripts/cbh_nc_indivual.py b/test_data/generate/cbh_nc_indivual.py similarity index 100% rename from test_data/scripts/cbh_nc_indivual.py rename to test_data/generate/cbh_nc_indivual.py diff --git a/test_data/scripts/cbh_to_nc.py b/test_data/generate/cbh_to_nc.py similarity index 100% rename from test_data/scripts/cbh_to_nc.py rename to test_data/generate/cbh_to_nc.py diff --git a/test_data/scripts/conftest.py b/test_data/generate/conftest.py similarity index 86% rename from test_data/scripts/conftest.py rename to test_data/generate/conftest.py index 2ad89d89..31468dd7 100644 --- a/test_data/scripts/conftest.py +++ b/test_data/generate/conftest.py @@ -32,17 +32,16 @@ def pytest_addoption(parser): @pytest.fixture() def exe(): - exe_name = "prms" platform = sys.platform.lower() if platform == "win32": - exe_name += "_win.exe" + exe_name = "prms_win_gfort_dbl_prec.exe" elif platform == "darwin": if processor() == "arm": - exe_name += "_mac_m1_intel" + exe_name = "prms_mac_m1_ifort_dbl_prec" else: - exe_name += "_mac" + exe_name = "prms_mac_intel_gfort_dbl_prec" elif platform == "linux": - exe_name += "_linux" + exe_name = "prms_linux_gfort_dbl_prec" exe_pth = pl.Path(f"../../bin/{exe_name}").resolve() return exe_pth @@ -161,3 +160,16 @@ def pytest_generate_tests(metafunc): ] ids = [pl.Path(kk).name for kk in simulations.keys()] metafunc.parametrize("simulation", sims, ids=ids, scope="session") + + if "final_file" in metafunc.fixturenames: + final_var_names = ["through_rain", "seg_lateral_inflow"] + final_files = [ + pl.Path(kk) / var + for kk in simulations.keys() + for var in final_var_names + ] + ids = [ff.parent.name + ":" + ff.name for ff in final_files] + # these are not really file names they are domain/key_var + metafunc.parametrize( + "final_file", final_files, ids=ids, scope="session" + ) diff --git a/test_data/scripts/conus_2yr_subset.sh b/test_data/generate/conus_2yr_subset.sh similarity index 100% rename from test_data/scripts/conus_2yr_subset.sh rename to test_data/generate/conus_2yr_subset.sh diff --git a/test_data/generate/convert_prms_output_to_nc.py b/test_data/generate/convert_prms_output_to_nc.py new file mode 100644 index 00000000..af560a6d --- /dev/null +++ b/test_data/generate/convert_prms_output_to_nc.py @@ -0,0 +1,58 @@ +from pathlib import Path +from filelock import FileLock + +import pytest + +from prms_convert_to_netcdf import convert_csv_to_nc, convert_soltab_to_nc +from prms_diagnostic_variables import ( + diagnose_simple_vars_to_nc, + diagnose_final_vars_to_nc, +) + + +@pytest.fixture +def netcdf_file(csv_file): + """Convert CSV files from model output to NetCDF""" + + var_name = csv_file.stem + data_dir = csv_file.parent + convert_csv_to_nc(var_name, data_dir) + + diagnose_simple_vars_to_nc(var_name, data_dir) + + return + + +def make_netcdf_files(netcdf_file): + print(f"Created NetCDF from CSV: {netcdf_file}") + + +@pytest.fixture() +def soltab_netcdf_file(tmp_path_factory, soltab_file): + """Convert soltab files to NetCDF, one file for each variable""" + domain_dir = soltab_file.parent + output_dir = domain_dir / "output" + convert_soltab_to_nc(soltab_file, output_dir, domain_dir) + + +def make_soltab_netcdf_files(soltab_netcdf_file): + print(f"Creating NetCDF files for soltab file {soltab_netcdf_file}") + + +@pytest.fixture(scope="session") +def final_netcdf_file(tmp_path_factory, final_file): + """Create NetCDF files that depend on multiple other NetCDFs""" + + domain_dir = final_file.parent + var_name = final_file.name + data_dir = domain_dir / "output" + + root_tmpdir = tmp_path_factory.getbasetemp().parent + with FileLock(root_tmpdir / "final_nc.lock"): + yield # do this in session cleanup + + diagnose_final_vars_to_nc(var_name, data_dir, domain_dir) + + +def make_final_netcdf_files(final_netcdf_file): + print(f"Creating final NetCDF file {final_netcdf_file}") diff --git a/test_data/scripts/drb_2yr_check.sh b/test_data/generate/drb_2yr_check.sh similarity index 100% rename from test_data/scripts/drb_2yr_check.sh rename to test_data/generate/drb_2yr_check.sh diff --git a/test_data/scripts/drb_2yr_subset.sh b/test_data/generate/drb_2yr_subset.sh similarity index 100% rename from test_data/scripts/drb_2yr_subset.sh rename to test_data/generate/drb_2yr_subset.sh diff --git a/test_data/generate/prms_convert_to_netcdf.py b/test_data/generate/prms_convert_to_netcdf.py new file mode 100644 index 00000000..e3c526c3 --- /dev/null +++ b/test_data/generate/prms_convert_to_netcdf.py @@ -0,0 +1,54 @@ +import pathlib as pl + +import pywatershed as pws +from pywatershed import CsvFile, Soltab + +"""This module is intended to aid/make consistent conversions of PRMS files""" + + +def convert_csv_to_nc( + var_name: str, data_dir: pl.Path, output_dir: pl.Path = None +): + """Convert PRMS CSV files to netcdf. + + Args: + var_name: str name of the variable to create + data_dir: where the csv file is found and the netcdf file will be + written (could add argument to output to a differnt dir) + """ + if output_dir is None: + output_dir = data_dir + csv_path = data_dir / f"{var_name}.csv" + nc_path = output_dir / f"{var_name}.nc" + CsvFile(csv_path).to_netcdf(nc_path) + assert nc_path.exists() + + +def convert_soltab_to_nc( + soltab_file: pl.Path, output_dir: pl.Path, domain_dir: pl.Path +): + """Convert soltab files to NetCDF, one file for each variable + + The inputs files are "soltab_debug" specifically written by the PRMS5.2.1 + in the pywatershed repository. + + Args: + soltab_file: the pl.Path of the soltab_debug file. + output_dir: pl.Path where the netcdf file output will be written + domain_dir: defaults to the parent dir of soltab_file, the pl.Path + where domain files (parameter & control) for the domain can be + found + """ + if domain_dir is None: + domain_dir = soltab_file.parent + + params = pws.parameters.PrmsParameters.load(domain_dir / "myparam.param") + nhm_ids = params.parameters["nhm_id"] + + soltab = Soltab(soltab_file, nhm_ids=nhm_ids) + soltab.to_netcdf(output_dir=output_dir) + print(f"Created NetCDF files from soltab file {soltab_file}:") + + for var in soltab.variables: + nc_path = output_dir / f"{var}.nc" + assert nc_path.exists() diff --git a/test_data/generate/prms_diagnostic_variables.py b/test_data/generate/prms_diagnostic_variables.py new file mode 100644 index 00000000..12610eb3 --- /dev/null +++ b/test_data/generate/prms_diagnostic_variables.py @@ -0,0 +1,332 @@ +import pathlib as pl + +import pywatershed as pws +from pywatershed.constants import epsilon32, nan, zero + +import numpy as np +import xarray as xr + +"""This module is for generating PRMS diagnostic variables""" + +# For generating timeseries of previous states +previous_vars = { + "gwres_stor": pws.PRMSGroundwater, + "dprst_stor_hru": pws.PRMSRunoff, + "hru_impervstor": pws.PRMSRunoff, + "pref_flow_stor": pws.PRMSSoilzone, + "slow_stor": pws.PRMSSoilzone, + "soil_lower": pws.PRMSSoilzone, + "soil_moist": pws.PRMSSoilzone, + "soil_rechr": pws.PRMSSoilzone, + "ssres_stor": pws.PRMSSoilzone, + "freeh2o": pws.PRMSSnow, + "pk_ice": pws.PRMSSnow, +} + +prev_rename = { + "gwres_stor": "gwres_stor_old", +} + +change_rename = {} + + +def diagnose_simple_vars_to_nc( + var_name: str, + data_dir: pl.Path, + domain_dir: pl.Path = None, + output_dir: pl.Path = None, +): + """Diagnose variables from PRMS output with only dependencies on var_name + + These variables are "simple" because they do not have any dependencies + besides the "var_name" variable which is a PRMS output variable. + + This is not to say that PRMS dosent output diagnostic variables, it + just dosent output these natively. + + Args: + var_name: str name of the variable to create + data_dir: where the netcdf file will be written, also where to look + for necessary inputs to create the diagnostic variable. + domain_dir: defaults to the parent dir of output_dir, this is where + domain files (parameter & control) for the domain can be found + output_dir: defaults to data_dir unless otherwise specified + + """ + + if domain_dir is None: + domain_dir = data_dir.parent + + if output_dir is None: + output_dir = data_dir + + nc_path = data_dir / f"{var_name}.nc" + + if var_name in previous_vars.keys(): + # the _prev is the desired suffix but PRMS legacy is inconsistent + current = xr.open_dataarray(nc_path) + prev = current * nan + prev[:] = np.roll(current.values, 1, axis=0) + + # the final value (-1) was wrapped to the zeroth position + # get the initial conditions for the first time by initializing the + # model This works based on the control file, so could handle restart. + control = pws.Control.load(domain_dir / "control.test") + control.options = control.options | { + "input_dir": domain_dir / "output", + } + params = pws.parameters.PrmsParameters.load( + domain_dir / "myparam.param" + ) + proc_class = previous_vars[var_name] + inputs = {} + for input_name in proc_class.get_inputs(): + # taken from process._initialize_var + meta = pws.meta.find_variables([input_name]) + dims = [params.dims[vv] for vv in meta[input_name]["dims"]] + init_type = meta[input_name]["type"] + + if len(dims) == 1: + inputs[input_name] = np.full(dims, nan, dtype=init_type) + else: + inputs[input_name] = pws.base.timeseries.TimeseriesArray( + var_name=input_name, + control=control, + array=np.full(dims, nan, dtype=init_type), + time=np.arange( + np.datetime64("2017-01-01"), + np.datetime64("2017-01-03"), + ), + ) + + proc = proc_class( + control=control, + discretization=None, + parameters=params, + **inputs, + ) + # these are the init_values + prev[0, :] = proc[var_name] + del proc + + change = current - prev + + # write the previous file + if var_name in prev_rename.keys(): + nc_name = f"{prev_rename[var_name]}.nc" + else: + nc_name = f"{var_name}_prev.nc" + out_nc_path = output_dir / nc_name + prev.rename(out_nc_path.stem).to_netcdf(out_nc_path) + assert nc_path.exists() + + # write the change file + if var_name in change_rename.keys(): + nc_name = f"{change_rename[var_name]}.nc" + else: + nc_name = f"{var_name}_change.nc" + out_nc_path = output_dir / nc_name + change.rename(out_nc_path.stem).to_netcdf(out_nc_path) + assert nc_path.exists() + return nc_path + + if var_name in ["sroff", "ssres_flow", "gwres_flow"]: + params = pws.parameters.PrmsParameters.load( + domain_dir / "myparam.param" + ) + ds = xr.open_dataset(nc_path) + ds = ds.rename({var_name: f"{var_name}_vol"}) + ds = ds * params.data_vars["hru_in_to_cf"] + ds.to_netcdf(output_dir / f"{var_name}_vol.nc") + ds.close() + + if var_name == "infil": + params = pws.parameters.PrmsParameters.load( + domain_dir / "myparam.param" + ).parameters + imperv_frac = params["hru_percent_imperv"] + dprst_frac = params["dprst_frac"] + perv_frac = 1.0 - imperv_frac - dprst_frac + ds = xr.open_dataset(nc_path.with_suffix(".nc")) + ds = ds.rename(infil="infil_hru") + ds["infil_hru"] = ds["infil_hru"] * perv_frac + ds.to_netcdf(output_dir / "infil_hru.nc") + ds.close() + + +def diagnose_final_vars_to_nc( + var_name: str, + data_dir: pl.Path, + domain_dir: pl.Path = None, + output_dir: pl.Path = None, +): + """Diagnose variables from multiple PRMS outputs + + In this case "var_name" is not a PRMS output variable, the variable is a + new variable which requires much/all of the PRMS output to be already + present in netcdf output format. This is why this is a final diagnostic + + Currently only: "through_rain" + + Args: + var_name: str name of the variable to create, not a PRMS variable. + output_dir: where the netcdf file will be written, also where to look + for necessary inputs to create the diagnostic variable. + domain_dir: defaults to the parent dir of output_dir, this is where + domain files (parameter & control) for the domain can be found + + """ + if domain_dir is None: + domain_dir = data_dir.parent + + if output_dir is None: + output_dir = data_dir + + if var_name == "through_rain": + out_file = output_dir / "through_rain.nc" + + data_vars = [ + "net_ppt", + "pptmix_nopack", + "snowmelt", + "pkwater_equiv", + "snow_evap", + "net_snow", + "net_rain", + ] + + data = {} + for vv in data_vars: + data_file = data_dir / f"{vv}.nc" + data[vv] = xr.open_dataarray(data_file) + + nearzero = 1.0e-6 + + cond1 = data["net_ppt"] > zero + cond2 = data["pptmix_nopack"] != 0 + cond3 = data["snowmelt"] < nearzero + cond4 = data["pkwater_equiv"] < epsilon32 + cond5 = data["snow_evap"] < nearzero + cond6 = data["net_snow"] < nearzero + + through_rain = data["net_rain"] * zero + # these are in reverse order + through_rain[:] = np.where( + cond1 & cond3 & cond4 & cond6, data["net_rain"], zero + ) + through_rain[:] = np.where( + cond1 & cond3 & cond4 & cond5, data["net_ppt"], through_rain + ) + through_rain[:] = np.where( + cond1 & cond2, data["net_rain"], through_rain + ) + + through_rain.to_dataset(name="through_rain").to_netcdf(out_file) + through_rain.close() + + for vv in data_vars: + data[vv].close() + + assert out_file.exists() + + if var_name in [ + "channel_sroff_vol", + "channel_ssres_flow_vol", + "channel_gwres_flow_vol", + "seg_lateral_inflow", + ]: + if var_name != "seg_lateral_inflow": + return + + data_vars = [ + "sroff_vol", + "ssres_flow_vol", + "gwres_flow_vol", + "seg_outflow", + "seg_inflow", + ] + data = {} + for vv in data_vars: + data_file = data_dir / f"{vv}.nc" + data[vv] = xr.open_dataarray(data_file) + + control = pws.Control.load(domain_dir / "control.test") + s_per_time = control.time_step_seconds + params = pws.parameters.PrmsParameters.load( + domain_dir / "myparam.param" + ) + + ntimes = control.n_times + nseg = params.dims["nsegment"] + nhru = params.dims["nhru"] + hru_segment = params.parameters["hru_segment"] - 1 + to_segment = (params.parameters["tosegment"] - 1).astype("int64") + + outflow_mask = np.full((len(to_segment)), False) + for iseg in range(nseg): + if to_segment[iseg] < 0: + outflow_mask[iseg] = True + continue + + seg_stor_change = ( + data["seg_inflow"] - data["seg_outflow"] + ) * s_per_time + + channel_outflow_vol = ( + np.where(outflow_mask, data["seg_outflow"], zero) + ) * s_per_time + + seg_lateral_inflow = np.zeros((ntimes, nseg)) + channel_sroff_vol = np.zeros((ntimes, nhru)) + channel_ssres_flow_vol = np.zeros((ntimes, nhru)) + channel_gwres_flow_vol = np.zeros((ntimes, nhru)) + + for ihru in range(nhru): + iseg = hru_segment[ihru] + if iseg < 0: + # This is bad, selective handling of fluxes is not cool, + # mass is being discarded in a way that has to be coordinated + # with other parts of the code. + # This code shuold be removed evenutally. + continue + + else: + channel_sroff_vol[:, ihru] = data["sroff_vol"].values[:, ihru] + channel_ssres_flow_vol[:, ihru] = data[ + "ssres_flow_vol" + ].values[:, ihru] + channel_gwres_flow_vol[:, ihru] = data[ + "gwres_flow_vol" + ].values[:, ihru] + + # cubicfeet to cfs + lateral_inflow = ( + channel_sroff_vol[:, ihru] + + channel_ssres_flow_vol[:, ihru] + + channel_gwres_flow_vol[:, ihru] + ) / (s_per_time) + + seg_lateral_inflow[:, iseg] += lateral_inflow + + for hvar in [ + "channel_sroff_vol", + "channel_ssres_flow_vol", + "channel_gwres_flow_vol", + "seg_lateral_inflow", + "seg_stor_change", + "channel_outflow_vol", + ]: + if hvar in [ + "seg_lateral_inflow", + "seg_stor_change", + "channel_outflow_vol", + ]: + dum = data["seg_outflow"].rename(hvar) * np.nan + else: + dum = data["sroff_vol"].rename(hvar) * np.nan + + dum[:, :] = locals()[hvar] + out_file = data_dir / f"{hvar}.nc" + dum.to_netcdf(out_file) + print(out_file) + assert out_file.exists() diff --git a/test_data/scripts/pytest.ini b/test_data/generate/pytest.ini similarity index 89% rename from test_data/scripts/pytest.ini rename to test_data/generate/pytest.ini index 51b9641a..b99a9aef 100644 --- a/test_data/scripts/pytest.ini +++ b/test_data/generate/pytest.ini @@ -3,6 +3,5 @@ addopts = --order-dependencies python_files = test_*.py python_functions = - create_* make_* test_* \ No newline at end of file diff --git a/test_data/generate/remove_prms_csvs.py b/test_data/generate/remove_prms_csvs.py new file mode 100644 index 00000000..86a139e9 --- /dev/null +++ b/test_data/generate/remove_prms_csvs.py @@ -0,0 +1,6 @@ +csvs_to_keep = ["hru_ppt", "intcp_stor", "potet", "gwres_stor"] + + +def test_csv_to_netcdf(csv_file): + if csv_file.with_suffix("").name not in (csvs_to_keep): + csv_file.unlink() diff --git a/test_data/scripts/test_run_domains.py b/test_data/generate/run_prms_domains.py similarity index 100% rename from test_data/scripts/test_run_domains.py rename to test_data/generate/run_prms_domains.py diff --git a/test_data/scripts/submit_domains_tg.sh b/test_data/generate/submit_domains_tg.sh similarity index 100% rename from test_data/scripts/submit_domains_tg.sh rename to test_data/generate/submit_domains_tg.sh diff --git a/test_data/hru_1/control.test b/test_data/hru_1/control.test index e692a8b6..74ea87b3 100644 --- a/test_data/hru_1/control.test +++ b/test_data/hru_1/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -76 +81 #### nhruOut_freq 1 @@ -251,12 +251,13 @@ nhruOut_format 5 #### nhruOutVar_names -76 +81 4 albedo cap_infil_tot contrib_fraction dprst_evap_hru +dprst_insroff_hru dprst_seep_hru dprst_stor_hru dprst_sroff_hru @@ -288,7 +289,11 @@ net_snow newsnow orad_hru perv_actet +pk_def +pk_den +pk_depth pk_ice +pk_temp pkwater_ante pkwater_equiv potet @@ -343,7 +348,7 @@ nsegmentOutON_OFF nsegmentOutVars 1 1 -9 +10 #### nsegmentOut_freq 1 @@ -351,8 +356,9 @@ nsegmentOut_freq 1 #### nsegmentOutVar_names -9 +10 4 +seg_inflow seg_outflow seg_gwflow seg_sroff diff --git a/test_data/scripts/test_nc_domains.py b/test_data/scripts/test_nc_domains.py deleted file mode 100644 index 38cf3686..00000000 --- a/test_data/scripts/test_nc_domains.py +++ /dev/null @@ -1,165 +0,0 @@ -from pathlib import Path -from filelock import FileLock - -import numpy as np - -from pywatershed import CsvFile, Soltab -from pywatershed.parameters import PrmsParameters -from pywatershed.constants import epsilon64, zero - -import pytest -import xarray as xr - - -# For generating timeseries of previous states -previous_vars = [ - "dprst_stor_hru", - "freeh2o", - "hru_impervstor", - "pk_ice", - "pref_flow_stor", - "slow_stor", - "soil_lower", - "soil_moist", - "soil_rechr", - "ssres_stor", -] - - -@pytest.fixture -def netcdf_file(csv_file) -> Path: - """Convert CSV files from model output to NetCDF""" - - nc_path = csv_file.with_suffix(".nc") - data_dir = csv_file.parent - domain_dir = data_dir.parent - CsvFile(csv_file).to_netcdf(nc_path) - assert nc_path.exists() - - if csv_file.stem in previous_vars: - nc_name = csv_file.stem + "_prev.nc" - nc_path = data_dir / nc_name - csv = CsvFile({nc_path.stem: csv_file}) - csv._get_data() # why so private? - orig_dates = csv._data["date"].copy() - csv._data = np.roll(csv._data, 1, axis=0) - csv._data["date"] = orig_dates - # Here we will eventually want to supply the desired initial conditions - for hh in range(len(csv._data[0])): - if hh == 0: - continue - csv._data[0][hh] = np.zeros([1])[0] - csv.to_netcdf(nc_path) - assert nc_path.exists() - - if nc_path.stem in ["sroff", "ssres_flow", "gwres_flow"]: - params = PrmsParameters.load(domain_dir / "myparam.param") - var = nc_path.stem - ds = xr.open_dataset(nc_path) - ds = ds.rename({var: f"{var}_vol"}) - ds = ds * params.data_vars["hru_in_to_cf"] - ds.to_netcdf(data_dir / f"{var}_vol.nc") - ds.close() - - if nc_path.stem == "infil": - params = PrmsParameters.load(domain_dir / "myparam.param").parameters - imperv_frac = params["hru_percent_imperv"] - dprst_frac = params["dprst_frac"] - perv_frac = 1.0 - imperv_frac - dprst_frac - ds = xr.open_dataset(nc_path.with_suffix(".nc")) - ds = ds.rename(infil="infil_hru") - ds["infil_hru"] = ds["infil_hru"] * perv_frac - ds.to_netcdf(data_dir / "infil_hru.nc") - ds.close() - - return nc_path - - -def make_netcdf_files(netcdf_file): - print(f"Created NetCDF from CSV: {netcdf_file}") - - -@pytest.fixture(scope="session") -def soltab_netcdf_file(tmp_path_factory, soltab_file) -> Path: - """Convert soltab files to NetCDF, one file for each variable""" - - # the nhm_ids are not available in the solta_debug file currently, so get - # them from the domain parameters - domain_dir = soltab_file.parent - output_dir = domain_dir / "output" - output_dir.mkdir(exist_ok=True) - root_tmpdir = tmp_path_factory.getbasetemp().parent - nc_paths = [] - - with FileLock(root_tmpdir / "soltab_nc.lock"): - yield # postpone the work until session cleanup - - # this is a hack that should probably rely on the yaml if/when this - # fails - params = PrmsParameters.load(domain_dir / "myparam.param") - nhm_ids = params.parameters["nhm_id"] - soltab = Soltab(soltab_file, nhm_ids=nhm_ids) - - already_created = (output_dir / "soltab_potsw.nc").is_file() - if not already_created: - soltab.to_netcdf(output_dir=output_dir) - print(f"Created NetCDF files from soltab file {soltab_file}:") - - for var in soltab.variables: - nc_path = output_dir / f"{var}.nc" - nc_paths.append(nc_path) - assert nc_path.exists() - print(f"\t{nc_path}") - - -def make_soltab_netcdf_files(soltab_netcdf_file): - print(f"Creating NetCDF files for soltab file {soltab_netcdf_file}") - - -@pytest.fixture(scope="session") -def final_netcdf_file(tmp_path_factory, simulation) -> Path: - """Create the final NetCDF file (through_rain.nc) from other NetCDFs""" - - root_tmpdir = tmp_path_factory.getbasetemp().parent - data_dir = Path(simulation["ws"]) / "output" - data_dir.mkdir(exist_ok=True) - nc_path = data_dir / "through_rain.nc" - - with FileLock(root_tmpdir / "final_nc.lock"): - yield # do this in session cleanup - - if nc_path.is_file(): - return - - data_vars = [ - "net_rain", - "pk_ice_prev", - "freeh2o_prev", - "newsnow", - "pptmix_nopack", - ] - - data = {} - for vv in data_vars: - data_file = data_dir / f"{vv}.nc" - result = xr.open_dataset(data_file)[vv] - data[vv] = result - - try: - wh_through = ( - ((data["pk_ice_prev"] + data["freeh2o_prev"]) <= epsilon64) - & ~(data["newsnow"] == 1) - ) | (data["pptmix_nopack"] == 1) - - through_rain = data["net_rain"].copy() - through_rain[:] = np.where(wh_through, data["net_rain"], zero) - through_rain.to_dataset(name="through_rain").to_netcdf(nc_path) - through_rain.close() - assert nc_path.exists() - finally: - for vv in data_vars: - data[vv].close() - - -def make_final_netcdf_files(final_netcdf_file): - print(f"Creating final NetCDF file {final_netcdf_file}") diff --git a/test_data/ucb_2yr/control.test b/test_data/ucb_2yr/control.test index df624841..effd9667 100644 --- a/test_data/ucb_2yr/control.test +++ b/test_data/ucb_2yr/control.test @@ -238,7 +238,7 @@ nhruOutON_OFF nhruOutVars 1 1 -78 +83 #### nhruOut_freq 1 @@ -251,13 +251,14 @@ nhruOut_format 4 #### nhruOutVar_names -78 +83 4 albedo cap_infil_tot contrib_fraction dprst_seep_hru dprst_evap_hru +dprst_insroff_hru dprst_sroff_hru dprst_stor_hru freeh2o @@ -288,7 +289,11 @@ net_snow newsnow orad_hru perv_actet +pk_def +pk_den +pk_depth pk_ice +pk_temp pkwater_ante pkwater_equiv potet @@ -345,7 +350,7 @@ nsegmentOutON_OFF nsegmentOutVars 1 1 -9 +10 #### nsegmentOut_freq 1 @@ -353,8 +358,9 @@ nsegmentOut_freq 1 #### nsegmentOutVar_names -9 +10 4 +seg_inflow seg_outflow seg_gwflow seg_sroff