From 90e6e219315eac423c00a4b4479905fc3f922c2c Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:08:54 +0100 Subject: [PATCH 01/12] Add pre-commit hooks. --- .flake8 | 2 ++ .pre-commit-config.yaml | 9 +++++++++ README.md | 22 ++++++++++++++++++++++ requirements-dev.txt | 1 + 4 files changed, 34 insertions(+) create mode 100644 .flake8 create mode 100644 .pre-commit-config.yaml create mode 100644 requirements-dev.txt diff --git a/.flake8 b/.flake8 new file mode 100644 index 00000000..aa079ec5 --- /dev/null +++ b/.flake8 @@ -0,0 +1,2 @@ +[flake8] +max-line-length=120 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..03424f0e --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,9 @@ +repos: +- repo: https://github.com/psf/black + rev: stable + hooks: + - id: black +- repo: git@github.com:PyCQA/flake8.git + rev: 7.3.0 + hooks: + - id: flake8 diff --git a/README.md b/README.md index 8845a352..fa71ea04 100644 --- a/README.md +++ b/README.md @@ -164,3 +164,25 @@ make: *** [Makefile:48: elmer_grid.o] Error 1 ``` *Solution:* `sudo apt-get install libproj-dev` + +# Developer notes + +## Code formatting + +Automated checks for PEP8 compiance are implemented following [^1] with some modifications. Please install pre-commit first: + +```sh +pipx install pre-commit +``` + +Then, install the pre-commit hooks by running the following command from this folder: + +```sh +pre-commit install +``` + +The pre-commit hooks will now run automatically on every commit. + +---- + +[^1]: https://ljvmiranda921.github.io/notebook/2018/06/21/precommits-using-black-and-flake8/ diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 00000000..dcee22a5 --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1 @@ +pre-commit~=4.0 From b792fec2b53a21793097829aaf8b6770a1da0913 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:14:01 +0100 Subject: [PATCH 02/12] Add workflow for running pre-commit. * Expected to fail --- .github/workflows/pre-commit.yml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 .github/workflows/pre-commit.yml diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml new file mode 100644 index 00000000..b049b81d --- /dev/null +++ b/.github/workflows/pre-commit.yml @@ -0,0 +1,18 @@ +name: Run pre-commit +on: + - push + - pull_request + +jobs: + pre-commit: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: '3.x' + - name: Install pre-commit + run: pip install pre-commit + - name: Run pre-commit + run: pre-commit run --all-files From f2042e0c3cdf6c12eb8d0ff765ee6a2284af7d61 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:20:23 +0100 Subject: [PATCH 03/12] Use HTTPS. --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 03424f0e..629c35a4 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,7 +3,7 @@ repos: rev: stable hooks: - id: black -- repo: git@github.com:PyCQA/flake8.git +- repo: https://github.com/PyCQA/flake8.git rev: 7.3.0 hooks: - id: flake8 From af08095bf0c9abe52021221c192b3c8db7b5292b Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:22:44 +0100 Subject: [PATCH 04/12] Only run flake8 in CI --- .github/workflows/pre-commit.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml index b049b81d..f580a789 100644 --- a/.github/workflows/pre-commit.yml +++ b/.github/workflows/pre-commit.yml @@ -15,4 +15,4 @@ jobs: - name: Install pre-commit run: pip install pre-commit - name: Run pre-commit - run: pre-commit run --all-files + run: pre-commit run flake8 --all-files From e9239e11dfd0c9a2c6dc1481f8a63c8341e304b3 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:24:48 +0100 Subject: [PATCH 05/12] Add note for bypassing. --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index fa71ea04..9e59e8b3 100644 --- a/README.md +++ b/README.md @@ -183,6 +183,8 @@ pre-commit install The pre-commit hooks will now run automatically on every commit. +Note: You can bypass this check with `--no-verify`. Please note that the CI will also check for PEP8 compliance and fail if the code is not PEP8 compliant. Therefore, it is recommended to use the pre-commit hooks locally before pushing code to this repository. + ---- [^1]: https://ljvmiranda921.github.io/notebook/2018/06/21/precommits-using-black-and-flake8/ From 2dd3bd4867e0aa808bb6bf13c782511eae66e881 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:42:27 +0100 Subject: [PATCH 06/12] Add line length. --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 pyproject.toml diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..55ec8d78 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[tool.black] +line-length = 120 From 31a8311bef0cd1539e751c3fd134f1388b2d0b0b Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 15:43:00 +0100 Subject: [PATCH 07/12] Run black and apply trivial changes. --- src/coupling.py | 204 +++++---- src/dummy.py | 160 ++++---- src/ebfm/FINAL_create_restart_file.py | 8 +- src/ebfm/INIT.py | 415 ++++++++++--------- src/ebfm/LOOP_EBM.py | 82 ++-- src/ebfm/LOOP_EBM_LHF.py | 19 +- src/ebfm/LOOP_EBM_LWin.py | 2 +- src/ebfm/LOOP_EBM_LWout.py | 2 +- src/ebfm/LOOP_EBM_SHF.py | 3 +- src/ebfm/LOOP_EBM_SWin.py | 41 +- src/ebfm/LOOP_EBM_SWout.py | 10 +- src/ebfm/LOOP_EBM_insolation.py | 72 ++-- src/ebfm/LOOP_SNOW.py | 567 +++++++++++++------------- src/ebfm/LOOP_climate_forcing.py | 96 +++-- src/ebfm/LOOP_general_functions.py | 9 +- src/ebfm/LOOP_mass_balance.py | 12 +- src/ebfm/LOOP_write_to_file.py | 23 +- src/ebfm/grid.py | 7 +- src/elmer/mesh.py | 48 ++- src/elmer/parser.py | 53 ++- src/reader.py | 98 +++-- src/utils.py | 21 +- 22 files changed, 1045 insertions(+), 907 deletions(-) diff --git a/src/coupling.py b/src/coupling.py index 09aa8bbc..a779a5a7 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -6,20 +6,22 @@ from typing import List, Dict from elmer.mesh import Mesh as Grid # for now use an alias + # from ebfm.geometry import Grid # TODO: consider introducing a new data structure native to EBFM? import logging + logger = logging.getLogger(__name__) -field_template = ''' +field_template = """ field {name}: - source: - component: {comp} - grid: {grid} - timestep: {timestep} - metadata: {metadata} -''' +""" def days_to_iso(days: float) -> str: @@ -30,10 +32,12 @@ def days_to_iso(days: float) -> str: @returns ISO 8601 formatted string representing the time step """ import pandas as pd + dt = pd.Timedelta(days=days) return dt.isoformat() -class Coupler(): + +class Coupler: component_name: str # name of this component in the coupler configuration couple_to_elmer_ice: bool # whether to couple with Elmer/Ice @@ -45,7 +49,7 @@ def has_coupling(self) -> bool: logger.debug(f"Component {self.component_name} has_coupling={return_value}") return return_value - def __init__(self, component_name: str, coupler_config: Path=None): + def __init__(self, component_name: str, coupler_config: Path = None): """ Create Coupler object @@ -77,13 +81,13 @@ def __init__(self, component_name: str): class YACCoupler(Coupler): - interface : yac.YAC = None - component : yac.Component = None - grid : yac.UnstructuredGrid = None - corner_points : yac.Points = None + interface: yac.YAC = None + component: yac.Component = None + grid: yac.UnstructuredGrid = None + corner_points: yac.Points = None # TODO use a single dictionary for both source and target fields; provide properties to create dicts containing source_fields and target_fields by comparing to exchange_type - source_fields : Dict[str, yac.Field] = {} - target_fields : Dict[str, yac.Field] = {} + source_fields: Dict[str, yac.Field] = {} + target_fields: Dict[str, yac.Field] = {} def __init__(self, component_name: str, coupler_config: Path): logger.debug(f"YAC version is {yac.version()}") @@ -142,25 +146,37 @@ def construct_coupling(self, time: Dict[str, float]): def construct_coupling_elmer_ice(self, time: Dict[str, float]): assert self.couple_to_elmer_ice, "Cannot construct coupling to Elmer/Ice if 'couple_to_elmer_ice' is False." - FieldDefinition = namedtuple('FieldDefinition', ['exchange_type', 'name', 'metadata']) + FieldDefinition = namedtuple("FieldDefinition", ["exchange_type", "name", "metadata"]) # TODO: Get hard-coded data below from dummies/EBFM/ebfm-config.yaml field_definitions = [ - FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="T_ice", metadata="Near surface temperature at Ice surface (in K)"), + FieldDefinition( + exchange_type=yac.ExchangeType.SOURCE, + name="T_ice", + metadata="Near surface temperature at Ice surface (in K)", + ), FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="smb", metadata="??? (in ???)"), FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="runoff", metadata="Runoff (in ???)"), FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name="h", metadata="Surface height (in m)"), ] - timestep_value = days_to_iso(time['dt']) + timestep_value = days_to_iso(time["dt"]) collection_size = 1 # TODO: Dummy value for now; make configurable if needed - Timestep = namedtuple('Timestep', ['value', 'format']) - timestep = Timestep(value=timestep_value, - format=yac.TimeUnit.ISO_FORMAT) + Timestep = namedtuple("Timestep", ["value", "format"]) + timestep = Timestep(value=timestep_value, format=yac.TimeUnit.ISO_FORMAT) for field_definition in field_definitions: - field = yac.Field.create(field_definition.name, self.component, self.corner_points, collection_size, timestep.value, timestep.format) - self.interface.def_field_metadata(field.component_name, field.grid_name, field_definition.name, field_definition.metadata.encode("utf-8")); + field = yac.Field.create( + field_definition.name, + self.component, + self.corner_points, + collection_size, + timestep.value, + timestep.format, + ) + self.interface.def_field_metadata( + field.component_name, field.grid_name, field_definition.name, field_definition.metadata.encode("utf-8") + ) if field_definition.exchange_type == yac.ExchangeType.SOURCE: self.source_fields[field_definition.name] = field @@ -171,31 +187,57 @@ def construct_coupling_elmer_ice(self, time: Dict[str, float]): def construct_coupling_icon_atmo(self, time: Dict[str, float]): assert self.couple_to_icon_atmo, "Cannot construct coupling to ICON if 'couple_to_icon_atmo' is False." - FieldDefinition = namedtuple('FieldDefinition', ['exchange_type', 'name', 'metadata']) + FieldDefinition = namedtuple("FieldDefinition", ["exchange_type", "name", "metadata"]) # TODO: Get hard-coded data below from dummies/EBFM/ebfm-config.yaml field_definitions = [ # FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="albedo", metadata="Albedo of the ice surface (in ???)"), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='pr', metadata='Precipitation rate (in kg m-2 s-1)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='pr_snow', metadata='Precipitation rate of snow (in kg m-2 s-1)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='rsds', metadata='Downward shortwave radiation flux (in W m-2)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='rlds', metadata='Downward longwave radiation flux (in W m-2)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='sfcwind', metadata='Wind speed at surface (in m s-1)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='clt', metadata='Cloud cover (in fraction)'), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='tas', metadata='Temperature at surface (in K)'), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, name="pr", metadata="Precipitation rate (in kg m-2 s-1)" + ), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, + name="pr_snow", + metadata="Precipitation rate of snow (in kg m-2 s-1)", + ), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, + name="rsds", + metadata="Downward shortwave radiation flux (in W m-2)", + ), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, + name="rlds", + metadata="Downward longwave radiation flux (in W m-2)", + ), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, name="sfcwind", metadata="Wind speed at surface (in m s-1)" + ), + FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name="clt", metadata="Cloud cover (in fraction)"), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, name="tas", metadata="Temperature at surface (in K)" + ), # FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='huss', metadata='Specific humidity at surface (in kg kg-1)'), # FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='sfcPressure', metadata='Surface pressure (in Pa)') ] - timestep_value = days_to_iso(time['dt']) + timestep_value = days_to_iso(time["dt"]) collection_size = 1 # TODO: Dummy value for now; make configurable if needed - Timestep = namedtuple('Timestep', ['value', 'format']) - timestep = Timestep(value=timestep_value, - format=yac.TimeUnit.ISO_FORMAT) + Timestep = namedtuple("Timestep", ["value", "format"]) + timestep = Timestep(value=timestep_value, format=yac.TimeUnit.ISO_FORMAT) for field_definition in field_definitions: - field = yac.Field.create(field_definition.name, self.component, self.corner_points, collection_size, timestep.value, timestep.format) - self.interface.def_field_metadata(field.component_name, field.grid_name, field_definition.name, field_definition.metadata.encode("utf-8")); + field = yac.Field.create( + field_definition.name, + self.component, + self.corner_points, + collection_size, + timestep.value, + timestep.format, + ) + self.interface.def_field_metadata( + field.component_name, field.grid_name, field_definition.name, field_definition.metadata.encode("utf-8") + ) if field_definition.exchange_type == yac.ExchangeType.SOURCE: self.source_fields[field_definition.name] = field @@ -208,18 +250,25 @@ def construct_coupling_post_sync(self): for field in self.target_fields.values(): is_defined = self.interface.get_field_is_defined(field.component_name, field.grid_name, field.name) - assert(is_defined), f"Field '{field.name}' is not defined in YAC for component '{field.component_name}' and grid '{field.grid_name}'." + assert ( + is_defined + ), f"Field '{field.name}' is not defined in YAC for component '{field.component_name}' and grid '{field.grid_name}'." field_role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) - assert(field_role == yac.ExchangeType.TARGET), f"Field '{field.name}' is no TARGET as expected but a {field_role}." + assert ( + field_role == yac.ExchangeType.TARGET + ), f"Field '{field.name}' is no TARGET as expected but a {field_role}." field_info = self.field_information(field) logger.info(field_info) - def field_information(self, field: yac.Field) -> str: - src_comp, src_grid, src_field = self.interface.get_field_source(field.component_name, field.grid_name, field.name) + src_comp, src_grid, src_field = self.interface.get_field_source( + field.component_name, field.grid_name, field.name + ) src_field_timestep = self.interface.get_field_timestep(src_comp, src_grid, src_field) src_field_metadata = self.interface.get_field_metadata(src_comp, src_grid, src_field) - return field_template.format(name=field.name, comp=src_comp, grid=src_grid, timestep=src_field_timestep, metadata=src_field_metadata) + return field_template.format( + name=field.name, comp=src_comp, grid=src_grid, timestep=src_field_timestep, metadata=src_field_metadata + ) def _put(self, field_name: str, field_data: np.array): field: yac.Field = self.source_fields[field_name] @@ -239,40 +288,46 @@ def exchange_icon_atmo(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra """ assert self.couple_to_icon_atmo, "Cannot exchange data with ICON atmosphere if 'couple_to_icon_atmo' is False." - icon_fields = set([ - # 'albedo', - 'pr', - 'pr_snow', - 'rsds', - 'rlds', - 'clt', - 'sfcwind', - 'tas', - # 'huss', - # 'sfcPressure' - ]) + icon_fields = set( + [ + # 'albedo', + "pr", + "pr_snow", + "rsds", + "rlds", + "clt", + "sfcwind", + "tas", + # 'huss', + # 'sfcPressure' + ] + ) for field_name, field in self.source_fields.items(): logger.debug(f"Checking field {field_name} for ICON atmosphere exchange...") if field_name in icon_fields: role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) - assert role == yac.ExchangeType.SOURCE, f"Field '{field_name}' is not a target field for ICON exchange, but has role '{role}'." - logger.debug(f'Sending field {field_name} to ICON atmosphere...') + assert ( + role == yac.ExchangeType.SOURCE + ), f"Field '{field_name}' is not a target field for ICON exchange, but has role '{role}'." + logger.debug(f"Sending field {field_name} to ICON atmosphere...") self._put(field_name, put_data[field_name]) else: - logger.debug(f'Skipping field {field_name} as it is not part of ICON atmosphere exchange.') + logger.debug(f"Skipping field {field_name} as it is not part of ICON atmosphere exchange.") received_data = {} for field_name, field in self.target_fields.items(): logger.debug(f"Checking field {field_name} for ICON atmosphere exchange...") if field_name in icon_fields: role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) - assert role == yac.ExchangeType.TARGET, f"Field '{field_name}' is not a target field for ICON exchange, but has role '{role}'." - logger.debug(f'Receiving field {field_name} from ICON atmosphere...') + assert ( + role == yac.ExchangeType.TARGET + ), f"Field '{field_name}' is not a target field for ICON exchange, but has role '{role}'." + logger.debug(f"Receiving field {field_name} from ICON atmosphere...") # Only get the first element, since we only have collection size of 1 received_data[field_name] = self._get(field_name)[0] else: - logger.debug(f'Skipping field {field_name} as it is not part of ICON atmosphere exchange.') + logger.debug(f"Skipping field {field_name} as it is not part of ICON atmosphere exchange.") return received_data @@ -287,38 +342,45 @@ def exchange_elmer_ice(self, put_data: Dict[str, np.array]) -> Dict[str, np.arra assert self.couple_to_elmer_ice, "Cannot exchange data with Elmer/Ice if 'couple_to_elmer_ice' is False." # TODO: make configurable if needed - elmer_fields = set([ - 'smb', - 'T_ice', - 'runoff', - 'h', - # 'dhdx', - # 'dhdy', - ]) + elmer_fields = set( + [ + "smb", + "T_ice", + "runoff", + "h", + # 'dhdx', + # 'dhdy', + ] + ) for field_name, field in self.source_fields.items(): if field_name in elmer_fields: role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) - assert role == yac.ExchangeType.SOURCE, f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." + assert ( + role == yac.ExchangeType.SOURCE + ), f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." self._put(field_name, put_data[field_name]) received_data = {} for field_name, field in self.target_fields.items(): if field_name in elmer_fields: role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) - assert role == yac.ExchangeType.TARGET, f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." + assert ( + role == yac.ExchangeType.TARGET + ), f"Field '{field_name}' is not a target field for Elmer/Ice exchange, but has role '{role}'." # Only get the first element, since we only have collection size of 1 received_data[field_name] = self._get(field_name)[0] return received_data def finalize(self): - """Finalize the coupling interface - """ + """Finalize the coupling interface""" del self.interface -def init(coupler_config: Path, ebfm_coupling_config: Path, couple_with_icon_atmo: bool, couple_with_elmer_ice: bool) -> Coupler: +def init( + coupler_config: Path, ebfm_coupling_config: Path, couple_with_icon_atmo: bool, couple_with_elmer_ice: bool +) -> Coupler: """Create interface to the coupler and register component @param[in] path to global Coupler configuration file @@ -331,7 +393,7 @@ def init(coupler_config: Path, ebfm_coupling_config: Path, couple_with_icon_atmo # TODO: use ebfm_coupling_config for EBFM specific configuration? - component_name = 'ebfm' # TODO: get from ebfm_coupling_config? + component_name = "ebfm" # TODO: get from ebfm_coupling_config? coupler = YACCoupler(component_name, coupler_config) coupler.couple_to_icon_atmo = couple_with_icon_atmo @@ -341,6 +403,7 @@ def init(coupler_config: Path, ebfm_coupling_config: Path, couple_with_icon_atmo return coupler + def setup(coupler: Coupler, grid: Grid, time: Dict[str, float]): """Performs initialization operations after init and before enterint the time loop @@ -349,11 +412,12 @@ def setup(coupler: Coupler, grid: Grid, time: Dict[str, float]): @param[in] grid Grid used by EBFM where coupling happens @param[in] time dictionary with time parameters, e.g. {'tn': 12, 'dt': 0.125} """ - grid_name = "ebfm_grid" # TODO: get from ebfm_coupling_config? + grid_name = "ebfm_grid" # TODO: get from ebfm_coupling_config? coupler.add_grid(grid_name, grid) coupler.add_couples(time) + def finalize(coupler: Coupler): logger.debug("Finalizing coupling...") coupler.finalize() diff --git a/src/dummy.py b/src/dummy.py index d316589e..8e635ed1 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -14,7 +14,7 @@ log_levels = { - 'file':logging.DEBUG, # log level for logging to file + "file": logging.DEBUG, # log level for logging to file 0: logging.INFO, # log level for rank 0 # 1: logging.DEBUG, # to log other ranks to console define log level here } @@ -29,38 +29,44 @@ def main(): parser = argparse.ArgumentParser() - parser.add_argument("--couple-to-elmer-ice", - action="store_true", - help="Enable coupling with Elmer/Ice models via YAC") - - parser.add_argument("--couple-to-icon-atmo", - action="store_true", - help="Enable coupling with ICON via YAC") - - parser.add_argument("--coupler-config", - type=Path, - help="Path to the coupling configuration file (YAC coupler_config.yaml).") - - parser.add_argument("--elmer-mesh", - type=Path, - help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required.") - - parser.add_argument("--matlab-mesh", - type=Path, - help="Path to the MATLAB mesh file. Either --elmer-mesh or --matlab-mesh is required.") - - parser.add_argument("--netcdf-mesh", - type=Path, - help="Path to the NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.") - - parser.add_argument("--is-partitioned-elmer-mesh", - action="store_true", - help="Indicate if the provided Elmer mesh is partitioned for parallel runs.") - - parser.add_argument("--use-part", - type=int, - default=MPI.COMM_WORLD.rank + 1, - help="If using a partitioned Elmer mesh, allows to specify which partition ID to use for this run. If not provided, the MPI rank + 1 will be used as partition ID.") + parser.add_argument( + "--couple-to-elmer-ice", action="store_true", help="Enable coupling with Elmer/Ice models via YAC" + ) + + parser.add_argument("--couple-to-icon-atmo", action="store_true", help="Enable coupling with ICON via YAC") + + parser.add_argument( + "--coupler-config", type=Path, help="Path to the coupling configuration file (YAC coupler_config.yaml)." + ) + + parser.add_argument( + "--elmer-mesh", type=Path, help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required." + ) + + parser.add_argument( + "--matlab-mesh", + type=Path, + help="Path to the MATLAB mesh file. Either --elmer-mesh or --matlab-mesh is required.", + ) + + parser.add_argument( + "--netcdf-mesh", + type=Path, + help="Path to the NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.", + ) + + parser.add_argument( + "--is-partitioned-elmer-mesh", + action="store_true", + help="Indicate if the provided Elmer mesh is partitioned for parallel runs.", + ) + + parser.add_argument( + "--use-part", + type=int, + default=MPI.COMM_WORLD.rank + 1, + help="If using a partitioned Elmer mesh, allows to specify which partition ID to use for this run. If not provided, the MPI rank + 1 will be used as partition ID.", + ) args = parser.parse_args() @@ -70,9 +76,9 @@ def main(): logger.debug(f" {arg}: {val}") # Model setup & initialization - grid, time2, io, phys = INIT.init_config(args) - C = INIT.init_constants() - grid = INIT.init_grid(grid, io, args) + grid, time2, io, phys = INIT.init_config(args) + C = INIT.init_constants() + grid = INIT.init_grid(grid, io, args) OUT, IN, OUTFILE = INIT.init_initial_conditions(C, grid, io, time2) @@ -81,17 +87,17 @@ def main(): # TODO consider introducing an ebfm_adapter_config.yaml coupler = coupling.init( coupler_config=args.coupler_config, - ebfm_coupling_config=Path('dummies') / 'EBFM' / 'ebfm-config.yaml', + ebfm_coupling_config=Path("dummies") / "EBFM" / "ebfm-config.yaml", couple_with_icon_atmo=args.couple_to_icon_atmo, couple_with_elmer_ice=args.couple_to_elmer_ice, ) coupling.setup(coupler, grid["mesh"], time2) else: - coupler = coupling.NoCoupler('ebfm') + coupler = coupling.NoCoupler("ebfm") # Time-loop - logger.info('Entering time loop...') - for t in range(1, time2['tn'] + 1): + logger.info("Entering time loop...") + for t in range(1, time2["tn"] + 1): # Print time to screen time2 = LOOP_general_functions.print_time(t, time2) @@ -100,27 +106,29 @@ def main(): # Read and prepare climate input if coupler and coupler.couple_to_icon_atmo: # Exchange data with ICON - logger.info('Data exchange with ICON') - logger.debug('Started...') + logger.info("Data exchange with ICON") + logger.debug("Started...") data_to_icon = { - 'albedo': OUT['albedo'], + "albedo": OUT["albedo"], } data_from_icon = coupler.exchange_icon_atmo(data_to_icon) - logger.debug('Done.') - logger.debug('Received the following data from ICON:', data_from_icon) - - IN['P'] = data_from_icon['pr']*time2['dt']*C['dayseconds']*1e-3 # convert units from kg m-2 s-1 to m w.e. - IN['snow'] = data_from_icon['pr_snow'] - IN['SWin'] = data_from_icon['rsds'] - IN['LWin'] = data_from_icon['rlds'] - IN['C'] = data_from_icon['clt'] - IN['WS'] = data_from_icon['sfcwind'] - IN['T'] = data_from_icon['tas'] - IN['rain'] = IN['P'] - IN['snow'] # TODO: make this more flexible and configurable - IN['q'][:] = 0 # TODO: Read q from ICON instead and convert to RH - IN['Pres'][:] = 101500 # TODO: Read Pres from ICON instead + logger.debug("Done.") + logger.debug("Received the following data from ICON:", data_from_icon) + + IN["P"] = ( + data_from_icon["pr"] * time2["dt"] * C["dayseconds"] * 1e-3 + ) # convert units from kg m-2 s-1 to m w.e. + IN["snow"] = data_from_icon["pr_snow"] + IN["SWin"] = data_from_icon["rsds"] + IN["LWin"] = data_from_icon["rlds"] + IN["C"] = data_from_icon["clt"] + IN["WS"] = data_from_icon["sfcwind"] + IN["T"] = data_from_icon["tas"] + IN["rain"] = IN["P"] - IN["snow"] # TODO: make this more flexible and configurable + IN["q"][:] = 0 # TODO: Read q from ICON instead and convert to RH + IN["Pres"][:] = 101500 # TODO: Read Pres from ICON instead IN, OUT = LOOP_climate_forcing.main(C, grid, IN, t, time2, OUT, coupler) @@ -128,53 +136,55 @@ def main(): OUT = LOOP_EBM.main(C, OUT, IN, time2, grid, coupler) # Run snow & firn model - OUT = LOOP_SNOW.main(C, OUT, IN, time2['dt'], grid, phys) + OUT = LOOP_SNOW.main(C, OUT, IN, time2["dt"], grid, phys) # Calculate surface mass balance OUT = LOOP_mass_balance.main(OUT, IN, C) if coupler.couple_to_elmer_ice: # Exchange data with Elmer - logger.info('Data exchange with Elmer/Ice') - logger.debug('Started...') + logger.info("Data exchange with Elmer/Ice") + logger.debug("Started...") data_to_elmer = { - 'smb': OUT['smb'], - 'T_ice': OUT['T_ice'], - 'runoff': OUT['runoff'], + "smb": OUT["smb"], + "T_ice": OUT["T_ice"], + "runoff": OUT["runoff"], } data_from_elmer = coupler.exchange_elmer_ice(data_to_elmer) - logger.debug('Done.') - logger.debug('Received the following data from Elmer/Ice:', data_from_elmer) + logger.debug("Done.") + logger.debug("Received the following data from Elmer/Ice:", data_from_elmer) - IN['h'] = data_from_elmer['h'] - grid['z'] = IN['h'][0].ravel() + IN["h"] = data_from_elmer["h"] + grid["z"] = IN["h"][0].ravel() # TODO add gradient field later # IN['dhdx'] = data_from_elmer('dhdx') # IN['dhdy'] = data_from_elmer('dhdy') # Write output to files (only in uncoupled run and for unpartitioned grid) - if not grid['is_partitioned'] and not coupler.has_coupling: - assert (grid['input_type'] is GridInputType.MATLAB), \ - "Output writing currently only implemented for MATLAB input grids." + if not grid["is_partitioned"] and not coupler.has_coupling: + assert ( + grid["input_type"] is GridInputType.MATLAB + ), "Output writing currently only implemented for MATLAB input grids." io, OUTFILE = LOOP_write_to_file.main(OUTFILE, io, OUT, grid, t, time2, C) pass - elif grid['is_partitioned'] or coupler.has_coupling: - logger.warning('Skipping writing output to file for coupled or partitioned runs.') + elif grid["is_partitioned"] or coupler.has_coupling: + logger.warning("Skipping writing output to file for coupled or partitioned runs.") else: - logger.error('Unhandled case in output writing.') + logger.error("Unhandled case in output writing.") raise Exception("Unhandled case in output writing.") # Write restart file - if not grid['is_partitioned'] and not coupler.has_coupling: + if not grid["is_partitioned"] and not coupler.has_coupling: FINAL_create_restart_file.main(OUT, io) - logger.info('Time loop completed.') + logger.info("Time loop completed.") if coupler.has_coupling: coupling.finalize(coupler) - logger.info('Closing down EBFM.') + logger.info("Closing down EBFM.") + # Entry point for script execution if __name__ == "__main__": diff --git a/src/ebfm/FINAL_create_restart_file.py b/src/ebfm/FINAL_create_restart_file.py index 9e75e67b..c571b1f9 100644 --- a/src/ebfm/FINAL_create_restart_file.py +++ b/src/ebfm/FINAL_create_restart_file.py @@ -2,6 +2,7 @@ from netCDF4 import Dataset, date2num import numpy as np + def main(OUT, io): def create_boot_file(): @@ -51,7 +52,9 @@ def create_boot_file(): return True - OUT['timelastsnow_netCDF'] = date2num(OUT['timelastsnow'], units="days since 1970-01-01 00:00:00", calendar="gregorian") + OUT["timelastsnow_netCDF"] = date2num( + OUT["timelastsnow"], units="days since 1970-01-01 00:00:00", calendar="gregorian" + ) OUT = { "subZ": OUT["subZ"], "subW": OUT["subW"], @@ -66,6 +69,5 @@ def create_boot_file(): "alb_snow": OUT["alb_snow"], } - # Create the boot file - create_boot_file() \ No newline at end of file + create_boot_file() diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index d5f8686f..aa7941e6 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -17,8 +17,10 @@ from mpi4py import MPI import logging + logger = logging.getLogger(__name__) + def init_config(args: Namespace): """ Set model parameters, specify grid parameters, model time period, I/O, and physics settings. @@ -33,36 +35,36 @@ def init_config(args: Namespace): # Time parameters # --------------------------------------------------------------------- time = {} - time['ts'] = datetime.strptime('1-Jan-1979 00:00', '%d-%b-%Y %H:%M') # Start date and time - time['te'] = datetime.strptime('2-Jan-1979 09:00', '%d-%b-%Y %H:%M') # End date and time - time['dt'] = 0.125 # Time step in days + time["ts"] = datetime.strptime("1-Jan-1979 00:00", "%d-%b-%Y %H:%M") # Start date and time + time["te"] = datetime.strptime("2-Jan-1979 09:00", "%d-%b-%Y %H:%M") # End date and time + time["dt"] = 0.125 # Time step in days # Calculate the number of time steps - time['tn'] = int(round((time['te'] - time['ts']).total_seconds() / 86400 / time['dt'])) + 1 - time['dT_UTC'] = 1 # Time difference relative to UTC in hours + time["tn"] = int(round((time["te"] - time["ts"]).total_seconds() / 86400 / time["dt"])) + 1 + time["dT_UTC"] = 1 # Time difference relative to UTC in hours # --------------------------------------------------------------------- # Grid parameters # --------------------------------------------------------------------- grid = {} - grid['utmzone'] = 33 # UTM zone - grid['max_subZ'] = 0.1 # Maximum first layer thickness (m) - grid['nl'] = 50 # Number of vertical layers - grid['doubledepth'] = True # Double vertical layer depth at specified layers (True/False) - grid['split'] = np.array([15, 25, 35]) # Vertical layer numbers at which layer depth doubles + grid["utmzone"] = 33 # UTM zone + grid["max_subZ"] = 0.1 # Maximum first layer thickness (m) + grid["nl"] = 50 # Number of vertical layers + grid["doubledepth"] = True # Double vertical layer depth at specified layers (True/False) + grid["split"] = np.array([15, 25, 35]) # Vertical layer numbers at which layer depth doubles # --------------------------------------------------------------------- # Model physics # --------------------------------------------------------------------- phys = {} - phys['percolation'] = 'normal' # Water percolation scheme + phys["percolation"] = "normal" # Water percolation scheme # Options: # - 'bucket': tipping-bucket method (all water added at the surface) # - 'normal': normally distributed deep percolation # - 'linear': linearly distributed deep percolation # - 'uniform': uniformly distributed deep percolation - phys['snow_compaction'] = 'firn+snow' # Snow and firn compaction scheme + phys["snow_compaction"] = "firn+snow" # Snow and firn compaction scheme # Options: # - 'firn_only': apply Ligtenberg et al. (2011) for all snow and firn layers # - 'firn+snow': apply Ligtenberg et al. (2011) for firn and Kampenhout et al. (2017) for seasonal snow @@ -72,20 +74,19 @@ def init_config(args: Namespace): # --------------------------------------------------------------------- io = {} - io['homedir'] = os.getcwd() # Home directory - io['outdir'] = os.path.join(io['homedir'], 'Output') # Output directory - io['rebootdir'] = os.path.join(io['homedir'], 'Reboot') # Restart file directory - io['readbootfile'] = False # REBOOT: read initial conditions from file (True/False) - io['writebootfile'] = True # REBOOT: write file for rebooting (True/False) - io['bootfilein'] = 'boot_final.nc' # REBOOT: bootfile to be read - io['bootfileout'] = 'boot_final.nc' # REBOOT: bootfile to be written - io['freqout'] = 8 # OUTPUT: frequency of storing output (every n-th time-step) - io['output_type'] = 2 # Set output file type: 1 = binary files, 2 = netCDF file + io["homedir"] = os.getcwd() # Home directory + io["outdir"] = os.path.join(io["homedir"], "Output") # Output directory + io["rebootdir"] = os.path.join(io["homedir"], "Reboot") # Restart file directory + io["readbootfile"] = False # REBOOT: read initial conditions from file (True/False) + io["writebootfile"] = True # REBOOT: write file for rebooting (True/False) + io["bootfilein"] = "boot_final.nc" # REBOOT: bootfile to be read + io["bootfileout"] = "boot_final.nc" # REBOOT: bootfile to be written + io["freqout"] = 8 # OUTPUT: frequency of storing output (every n-th time-step) + io["output_type"] = 2 # Set output file type: 1 = binary files, 2 = netCDF file # Ensure output and reboot directories exist - os.makedirs(io['outdir'], exist_ok=True) - os.makedirs(io['rebootdir'], exist_ok=True) - + os.makedirs(io["outdir"], exist_ok=True) + os.makedirs(io["rebootdir"], exist_ok=True) # Return the initialized parameters return grid, time, io, phys @@ -103,56 +104,56 @@ def init_constants(): # --------------------------------------------------------------------- # Energy balance model # --------------------------------------------------------------------- - C['alb_fresh'] = 0.83 # Albedo fresh snow (fraction) - C['alb_ice'] = 0.39 # Albedo ice (fraction) - C['alb_firn'] = 0.52 # Albedo firn (fraction) - C['tstar_wet'] = 15 # Albedo decay time-scale wet snow (days) - C['tstar_dry'] = 30 # Albedo decay time-scale dry snow (days) - C['tstar_K'] = 7 # Albedo decay time-scale coefficient - C['dstar'] = 7.0 # Albedo characteristic depth (mm w.e.) - C['b'] = 0.455 # Constant in LWin formulation - C['ecl'] = 0.960 # Constant in LWin formulation - C['p'] = 2 # Exponent in LWin formulation - C['boltz'] = 5.67e-8 # Stefan-Boltzmann constant (W m-2 K-4) - C['VP0'] = 610.5 # Vapor pressure at 0 degrees C - C['Cp'] = 1005.7 # Specific heat of dry air (J kg-1 K-1) - C['Cw'] = 4187.0 # Specific heat of water (J kg-1 K-1) - C['Ls'] = 2.83e6 # Latent heat of sublimation/riming (J kg-1) - C['Lm'] = 0.33e6 # Latent heat of melting/fusion (J kg-1) - C['Lv'] = 2.5e6 # Latent heat of evaporation/condensation (J kg-1) - C['Rv'] = 462.0 # Specific gas constant water vapor (J kg-1 K-1) - C['Rd'] = 287.0 # Specific gas constant dry air (J kg-1 K-1) - C['eps'] = 0.622 # C.Rd/C.Rv - C['dTacc'] = 0.01 # Threshold dT in solving the energy balance equation (K) - C['Pref'] = 1015e2 # Reference air pressure (Pa) - C['Pr'] = 5 # Prandtl number in SHF/LHF formulation - C['T0'] = 273.15 # Melting temperature of ice (K) - C['g'] = 9.81 # Gravitational acceleration (m s-2) - C['rd'] = 8.314 # Universal gas constant (J mol-1 K-1) - C['k_aer'] = 0.974 # Aerosol transmissivity exponent - C['k_turb'] = 0.0004 # Turbulent flux coefficient - C['turb'] = 0.0025 # Background turbulent exchange coefficient - C['rainsnowT'] = 273.75 # Temperature of snow to rain transition (K) - C['Pthres'] = 2.5e-8 # Threshold precipitation to reset time since last snowfall (m w.e. s-1) + C["alb_fresh"] = 0.83 # Albedo fresh snow (fraction) + C["alb_ice"] = 0.39 # Albedo ice (fraction) + C["alb_firn"] = 0.52 # Albedo firn (fraction) + C["tstar_wet"] = 15 # Albedo decay time-scale wet snow (days) + C["tstar_dry"] = 30 # Albedo decay time-scale dry snow (days) + C["tstar_K"] = 7 # Albedo decay time-scale coefficient + C["dstar"] = 7.0 # Albedo characteristic depth (mm w.e.) + C["b"] = 0.455 # Constant in LWin formulation + C["ecl"] = 0.960 # Constant in LWin formulation + C["p"] = 2 # Exponent in LWin formulation + C["boltz"] = 5.67e-8 # Stefan-Boltzmann constant (W m-2 K-4) + C["VP0"] = 610.5 # Vapor pressure at 0 degrees C + C["Cp"] = 1005.7 # Specific heat of dry air (J kg-1 K-1) + C["Cw"] = 4187.0 # Specific heat of water (J kg-1 K-1) + C["Ls"] = 2.83e6 # Latent heat of sublimation/riming (J kg-1) + C["Lm"] = 0.33e6 # Latent heat of melting/fusion (J kg-1) + C["Lv"] = 2.5e6 # Latent heat of evaporation/condensation (J kg-1) + C["Rv"] = 462.0 # Specific gas constant water vapor (J kg-1 K-1) + C["Rd"] = 287.0 # Specific gas constant dry air (J kg-1 K-1) + C["eps"] = 0.622 # C.Rd/C.Rv + C["dTacc"] = 0.01 # Threshold dT in solving the energy balance equation (K) + C["Pref"] = 1015e2 # Reference air pressure (Pa) + C["Pr"] = 5 # Prandtl number in SHF/LHF formulation + C["T0"] = 273.15 # Melting temperature of ice (K) + C["g"] = 9.81 # Gravitational acceleration (m s-2) + C["rd"] = 8.314 # Universal gas constant (J mol-1 K-1) + C["k_aer"] = 0.974 # Aerosol transmissivity exponent + C["k_turb"] = 0.0004 # Turbulent flux coefficient + C["turb"] = 0.0025 # Background turbulent exchange coefficient + C["rainsnowT"] = 273.75 # Temperature of snow to rain transition (K) + C["Pthres"] = 2.5e-8 # Threshold precipitation to reset time since last snowfall (m w.e. s-1) # --------------------------------------------------------------------- # Snow model # --------------------------------------------------------------------- - C['Dfreshsnow'] = 350.0 # Density of fresh snow (kg m-3) - C['Dice'] = 900.0 # Density of ice (kg m-3) - C['Dfirn'] = 500.0 # Density of firn (kg m-3) - C['Dwater'] = 1000.0 # Density of water (kg m-3) - C['Ec'] = 60000 # Gravitational densification factor - C['Eg'] = 42400 # Gravitational densification factor - C['Trunoff'] = 0.001 # Slush runoff time-scale (days) - C['geothermal_flux'] = 0.0 # Geothermal heat flux (W m-2) - C['perc_depth'] = 6.0 # Characteristic deep percolation depth (m) + C["Dfreshsnow"] = 350.0 # Density of fresh snow (kg m-3) + C["Dice"] = 900.0 # Density of ice (kg m-3) + C["Dfirn"] = 500.0 # Density of firn (kg m-3) + C["Dwater"] = 1000.0 # Density of water (kg m-3) + C["Ec"] = 60000 # Gravitational densification factor + C["Eg"] = 42400 # Gravitational densification factor + C["Trunoff"] = 0.001 # Slush runoff time-scale (days) + C["geothermal_flux"] = 0.0 # Geothermal heat flux (W m-2) + C["perc_depth"] = 6.0 # Characteristic deep percolation depth (m) # --------------------------------------------------------------------- # Other constants # --------------------------------------------------------------------- - C['yeardays'] = 365.242199 # Days in a year - C['dayseconds'] = 24 * 3600 # Seconds per day + C["yeardays"] = 365.242199 # Days in a year + C["dayseconds"] = 24 * 3600 # Seconds per day return C @@ -169,33 +170,35 @@ def init_grid(grid, io, args: Namespace): logger.error("--is-partitioned-elmer-mesh requires --elmer-mesh.") raise Exception("Invalid grid configuration.") - grid['is_partitioned'] = args.is_partitioned_elmer_mesh - if grid['is_partitioned']: - assert (args.netcdf_mesh), "--is-partitioned-elmer-mesh requires --netcdf-mesh. (Without --netcdf-mesh should also work but is untested.)" - logger.info('Using partitioned grid...') + grid["is_partitioned"] = args.is_partitioned_elmer_mesh + if grid["is_partitioned"]: + assert ( + args.netcdf_mesh + ), "--is-partitioned-elmer-mesh requires --netcdf-mesh. (Without --netcdf-mesh should also work but is untested.)" + logger.info("Using partitioned grid...") else: - logger.info('Using non-partitioned grid...') + logger.info("Using non-partitioned grid...") if args.matlab_mesh: - grid['input_type'] = GridInputType.MATLAB + grid["input_type"] = GridInputType.MATLAB elif args.netcdf_mesh and args.elmer_mesh: - grid['input_type'] = GridInputType.CUSTOM + grid["input_type"] = GridInputType.CUSTOM elif args.elmer_mesh: - grid['input_type'] = GridInputType.ELMER + grid["input_type"] = GridInputType.ELMER else: - logger.error(f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. Please please refer to the documentation for correct configuration.") + logger.error( + f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. Please please refer to the documentation for correct configuration." + ) raise Exception("Invalid grid configuration.") - if grid['input_type'] is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine - if grid['is_partitioned']: - mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh, - is_partitioned=True, - partition_id=args.use_part) + if grid["input_type"] is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine + if grid["is_partitioned"]: + mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh, is_partitioned=True, partition_id=args.use_part) else: mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh) - grid['x'], grid['y'] = mesh.x_vertices, mesh.y_vertices - grid['z'] = read_dem(args.netcdf_mesh, grid['x'], grid['y']) + grid["x"], grid["y"] = mesh.x_vertices, mesh.y_vertices + grid["z"] = read_dem(args.netcdf_mesh, grid["x"], grid["y"]) grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! @@ -203,122 +206,125 @@ def init_grid(grid, io, args: Namespace): grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! grid["slope_gamma"] = np.zeros_like(grid["x"]) # test values! grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier - grid["gpsum"] = np.sum(grid['mask'] == 1) # number of modelled grid cells + grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells grid["mesh"] = mesh # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy - elif grid['input_type'] is GridInputType.ELMER: # Read grid and elevations from Elmer - mesh: Mesh = read_elmer_mesh(args.elmer_mesh) # assuming mesh/MESH/mesh.nodes contains DEM data in the z component; see mesh/README.md for the required preprocessing steps. + elif grid["input_type"] is GridInputType.ELMER: # Read grid and elevations from Elmer + mesh: Mesh = read_elmer_mesh( + args.elmer_mesh + ) # assuming mesh/MESH/mesh.nodes contains DEM data in the z component; see mesh/README.md for the required preprocessing steps. grid["x"], grid["y"], grid["z"] = mesh.x_vertices, mesh.y_vertices, mesh.z_vertices - grid["z"] = np.random.uniform(0, 100, size=len(grid['x'])) # test values! - grid["slope_x"] = np.zeros_like(grid["x"]) # test values! - grid["slope_y"] = np.zeros_like(grid["x"]) # test values! - grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! - grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! - grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! + grid["z"] = np.random.uniform(0, 100, size=len(grid["x"])) # test values! + grid["slope_x"] = np.zeros_like(grid["x"]) # test values! + grid["slope_y"] = np.zeros_like(grid["x"]) # test values! + grid["lat"] = np.zeros_like(grid["x"]) + 75 # test values! + grid["lon"] = np.zeros_like(grid["x"]) + 320 # test values! + grid["slope_beta"] = np.zeros_like(grid["x"]) # test values! grid["slope_gamma"] = np.zeros_like(grid["x"]) # test values! - grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier - grid["gpsum"] = np.sum(grid['mask'] == 1) # number of modelled grid cells + grid["mask"] = np.ones_like(grid["x"]) # treats every grid cell as glacier + grid["gpsum"] = np.sum(grid["mask"] == 1) # number of modelled grid cells # TODO later add slope # grid["slope_x"], grid["slope_y"] = mesh.dzdy, mesh.dzdy grid["mesh"] = mesh - elif grid['input_type'] is GridInputType.MATLAB: # Read grid and elevations from example MATLAB file + elif grid["input_type"] is GridInputType.MATLAB: # Read grid and elevations from example MATLAB file # --------------------------------------------------------------------- # Read and process grid information # --------------------------------------------------------------------- # Read grid data input_data = read_MATLAB_grid(args.matlab_mesh) - grid['x_2D'] = input_data['x'][0][0] - grid['y_2D'] = input_data['y'][0][0] - grid['z_2D'] = input_data['z'][0][0] - grid['mask_2D'] = input_data['mask'][0][0] + grid["x_2D"] = input_data["x"][0][0] + grid["y_2D"] = input_data["y"][0][0] + grid["z_2D"] = input_data["z"][0][0] + grid["mask_2D"] = input_data["mask"][0][0] # Determine domain extent - grid['Lx'], grid['Ly'] = grid['x_2D'].shape + grid["Lx"], grid["Ly"] = grid["x_2D"].shape # Flip grid E-W or N-S when needed - fy, _ = np.gradient(grid['y_2D']) + fy, _ = np.gradient(grid["y_2D"]) if fy[0, 0] < 0: - grid['x_2D'] = np.flipud(grid['x_2D']) - grid['y_2D'] = np.flipud(grid['y_2D']) - grid['z_2D'] = np.flipud(grid['z_2D']) - grid['mask_2D'] = np.flipud(grid['mask_2D']) + grid["x_2D"] = np.flipud(grid["x_2D"]) + grid["y_2D"] = np.flipud(grid["y_2D"]) + grid["z_2D"] = np.flipud(grid["z_2D"]) + grid["mask_2D"] = np.flipud(grid["mask_2D"]) - _, fx = np.gradient(grid['x_2D']) + _, fx = np.gradient(grid["x_2D"]) if fx[0, 0] < 0: - grid['x_2D'] = np.fliplr(grid['x_2D']) - grid['y_2D'] = np.fliplr(grid['y_2D']) - grid['z_2D'] = np.fliplr(grid['z_2D']) - grid['mask_2D'] = np.fliplr(grid['mask_2D']) + grid["x_2D"] = np.fliplr(grid["x_2D"]) + grid["y_2D"] = np.fliplr(grid["y_2D"]) + grid["z_2D"] = np.fliplr(grid["z_2D"]) + grid["mask_2D"] = np.fliplr(grid["mask_2D"]) # Calculate grid spacing - grid['dx'] = grid['x_2D'][0][1] - grid['x_2D'][0][0] + grid["dx"] = grid["x_2D"][0][1] - grid["x_2D"][0][0] # Calculate number of modeled grid cells (gpsum) - grid['gpsum'] = np.sum(grid['mask_2D'] == 1) + grid["gpsum"] = np.sum(grid["mask_2D"] == 1) # Create 1-D mask - grid['mask'] = grid['mask_2D'][grid['mask_2D'] == 1] + grid["mask"] = grid["mask_2D"][grid["mask_2D"] == 1] # Calculate latitude & longitude fields (from the original UTM coordinates) - utmzone = grid['utmzone'] # Assume this is already part of the grid + utmzone = grid["utmzone"] # Assume this is already part of the grid utm_to_latlon = Transformer.from_crs(f"EPSG:{32600 + utmzone}", "EPSG:4326", always_xy=True) - x_coords = grid['x_2D'].ravel() - y_coords = grid['y_2D'].ravel() + x_coords = grid["x_2D"].ravel() + y_coords = grid["y_2D"].ravel() lon, lat = utm_to_latlon.transform(x_coords, y_coords) # Reshape to 2D arrays matching the input's shape - grid['lon_2D'] = lon.reshape(grid['x_2D'].shape) - grid['lat_2D'] = lat.reshape(grid['y_2D'].shape) + grid["lon_2D"] = lon.reshape(grid["x_2D"].shape) + grid["lat_2D"] = lat.reshape(grid["y_2D"].shape) # Store 1-D (vectorized) grid information - mask_flat = grid['mask_2D'].flatten() - grid['x'] = grid['x_2D'].flatten()[mask_flat == 1] - grid['y'] = grid['y_2D'].flatten()[mask_flat == 1] - grid['z'] = grid['z_2D'].flatten()[mask_flat == 1] - grid['ind'] = np.where(grid['mask_2D'].flatten() == 1) - grid['xind'], grid['yind'] = np.where(grid['mask_2D'] == 1) + mask_flat = grid["mask_2D"].flatten() + grid["x"] = grid["x_2D"].flatten()[mask_flat == 1] + grid["y"] = grid["y_2D"].flatten()[mask_flat == 1] + grid["z"] = grid["z_2D"].flatten()[mask_flat == 1] + grid["ind"] = np.where(grid["mask_2D"].flatten() == 1) + grid["xind"], grid["yind"] = np.where(grid["mask_2D"] == 1) # --------------------------------------------------------------------- # Grid slope and aspect # --------------------------------------------------------------------- # Calculate slope and aspect - gradN, gradE = np.gradient(grid['z_2D'], grid['y_2D'][:, 0], grid['x_2D'][0, :]) - slope_rad = np.arctan(np.sqrt(gradN ** 2 + gradE ** 2)) + gradN, gradE = np.gradient(grid["z_2D"], grid["y_2D"][:, 0], grid["x_2D"][0, :]) + slope_rad = np.arctan(np.sqrt(gradN**2 + gradE**2)) slope_deg = np.degrees(slope_rad) aspect = np.arctan2(-gradE, -gradN) * (180 / np.pi) aspect[aspect < 0] = 360 + aspect[aspect < 0] - grid['slope'] = np.tan(np.radians(slope_deg)) - grid['slope_x'] = gradE - grid['slope_y'] = gradN - grid['aspect'] = aspect - grid['slope_2D'] = grid['slope'] - grid['slope_beta_2D'] = np.arctan(grid['slope']) + grid["slope"] = np.tan(np.radians(slope_deg)) + grid["slope_x"] = gradE + grid["slope_y"] = gradN + grid["aspect"] = aspect + grid["slope_2D"] = grid["slope"] + grid["slope_beta_2D"] = np.arctan(grid["slope"]) # Convert 2-D fields to 1-D vectors - grid['slope'] = grid['slope'].flatten()[mask_flat == 1] - grid['slope_x'] = grid['slope_x'].flatten()[mask_flat == 1] - grid['slope_y'] = grid['slope_y'].flatten()[mask_flat == 1] - grid['aspect'] = grid['aspect'].flatten()[mask_flat == 1] - grid['lat'] = grid['lat_2D'].flatten()[mask_flat == 1] - grid['lon'] = grid['lon_2D'].flatten()[mask_flat == 1] + grid["slope"] = grid["slope"].flatten()[mask_flat == 1] + grid["slope_x"] = grid["slope_x"].flatten()[mask_flat == 1] + grid["slope_y"] = grid["slope_y"].flatten()[mask_flat == 1] + grid["aspect"] = grid["aspect"].flatten()[mask_flat == 1] + grid["lat"] = grid["lat_2D"].flatten()[mask_flat == 1] + grid["lon"] = grid["lon_2D"].flatten()[mask_flat == 1] # Calculate slope_beta and slope_gamma (defining the tilt and orientation of a sloping surface) - grid['slope_beta'] = np.arctan(grid['slope']) - grid['slope_gamma'] = np.zeros_like(grid['slope']) - grid['slope_gamma'][grid['slope_y'] > 0] = np.arctan( - -grid['slope_x'][grid['slope_y'] > 0] / grid['slope_y'][grid['slope_y'] > 0]) - cond = np.logical_and(grid['slope_y'] < 0, grid['slope_x'] > 0) - grid['slope_gamma'][cond] = -np.pi + np.arctan(-grid['slope_x'][cond] / grid['slope_y'][cond]) - cond = np.logical_and(grid['slope_y'] < 0, grid['slope_x'] < 0) - grid['slope_gamma'][cond] = np.pi + np.arctan(-grid['slope_x'][cond] / grid['slope_y'][cond]) - grid['slope_gamma'][(grid['slope_x'] == 0) & (grid['slope_y'] < 0)] = np.pi - grid['slope_gamma'][(grid['slope_x'] == 0) & (grid['slope_y'] == 0)] = 0 - grid['slope_gamma'][(grid['slope_x'] > 0) & (grid['slope_y'] == 0)] = np.pi / 2 - grid['slope_gamma'][(grid['slope_x'] < 0) & (grid['slope_y'] == 0)] = -np.pi / 2 - grid['slope_gamma'] = -grid['slope_gamma'] + grid["slope_beta"] = np.arctan(grid["slope"]) + grid["slope_gamma"] = np.zeros_like(grid["slope"]) + grid["slope_gamma"][grid["slope_y"] > 0] = np.arctan( + -grid["slope_x"][grid["slope_y"] > 0] / grid["slope_y"][grid["slope_y"] > 0] + ) + cond = np.logical_and(grid["slope_y"] < 0, grid["slope_x"] > 0) + grid["slope_gamma"][cond] = -np.pi + np.arctan(-grid["slope_x"][cond] / grid["slope_y"][cond]) + cond = np.logical_and(grid["slope_y"] < 0, grid["slope_x"] < 0) + grid["slope_gamma"][cond] = np.pi + np.arctan(-grid["slope_x"][cond] / grid["slope_y"][cond]) + grid["slope_gamma"][(grid["slope_x"] == 0) & (grid["slope_y"] < 0)] = np.pi + grid["slope_gamma"][(grid["slope_x"] == 0) & (grid["slope_y"] == 0)] = 0 + grid["slope_gamma"][(grid["slope_x"] > 0) & (grid["slope_y"] == 0)] = np.pi / 2 + grid["slope_gamma"][(grid["slope_x"] < 0) & (grid["slope_y"] == 0)] = -np.pi / 2 + grid["slope_gamma"] = -grid["slope_gamma"] return grid @@ -338,7 +344,7 @@ def read_MATLAB_grid(gridfile: Path): - 'mask': 2D NumPy array of glacier mask (0 = no glacier, 1 = glacier). """ - logger.info('EBFM: Reading grid data from MATLAB file...') + logger.info("EBFM: Reading grid data from MATLAB file...") input_data: dict[str, ndarray[tuple[int, ...], dtype[Any]]] = {} try: @@ -383,16 +389,16 @@ def init_initial_conditions(C, grid, io, time): OUT = {} # Dictionary to hold model output variables IN = {} # Dictionary to hold model input variables - gpsum = grid['gpsum'] - nl = grid['nl'] + gpsum = grid["gpsum"] + nl = grid["nl"] ########################################################## # Initialize conditions from restart file or set manually ########################################################## - if io.get('readbootfile', False): - logger.info('EBFM: Initialize from restart file...') + if io.get("readbootfile", False): + logger.info("EBFM: Initialize from restart file...") - reboot_dir = io['rebootdir'] + reboot_dir = io["rebootdir"] boot_filepath = f"{reboot_dir}/{io['bootfilein']}" # Open the NetCDF file @@ -407,66 +413,67 @@ def init_initial_conditions(C, grid, io, time): # Save the variable data to the dictionary OUT[var_name] = var_data - OUT['timelastsnow'] = num2date(OUT['timelastsnow_netCDF'], units="days since 1970-01-01 00:00:00", calendar="gregorian") + OUT["timelastsnow"] = num2date( + OUT["timelastsnow_netCDF"], units="days since 1970-01-01 00:00:00", calendar="gregorian" + ) else: - logger.info('EBFM: Initialize from manually set conditions...') - - OUT['Tsurf'] = np.full((gpsum,), 273.15) # Surface temperature (K) - OUT['subT'] = np.full((gpsum, nl), 265.0) # Vertical temperatures (K) - OUT['subW'] = np.zeros((gpsum, nl)) # Vertical irreducible water content (kg) - OUT['subS'] = np.zeros((gpsum, nl)) # Vertical slush water content (kg) - OUT['subD'] = np.full((gpsum, nl), C['Dice']) # Vertical densities (kg m-3) - OUT['subTmean'] = OUT['subT'] # Annual mean vertical layer temperature (K) - OUT['timelastsnow'] = np.full((gpsum,), time['ts']) # Timestep of last snowfall (days) - OUT['ys'] = np.full((gpsum,), 500.0) # Annual snowfall (mm water equivalent) - OUT['subZ'] = np.full((gpsum, nl), grid['max_subZ']) # Vertical layer depths (m) - OUT['alb_snow'] = np.full((gpsum,), C['alb_fresh']) # Snow albedo - OUT['snowmass'] = np.zeros((gpsum,)) # Snow mass (m water equivalent) - - if grid.get('doubledepth', False): # Sets layer thicknesses when 'double depth' is active - mask_indices = np.where(grid['mask'] == 1) - split = grid['split'] + logger.info("EBFM: Initialize from manually set conditions...") + + OUT["Tsurf"] = np.full((gpsum,), 273.15) # Surface temperature (K) + OUT["subT"] = np.full((gpsum, nl), 265.0) # Vertical temperatures (K) + OUT["subW"] = np.zeros((gpsum, nl)) # Vertical irreducible water content (kg) + OUT["subS"] = np.zeros((gpsum, nl)) # Vertical slush water content (kg) + OUT["subD"] = np.full((gpsum, nl), C["Dice"]) # Vertical densities (kg m-3) + OUT["subTmean"] = OUT["subT"] # Annual mean vertical layer temperature (K) + OUT["timelastsnow"] = np.full((gpsum,), time["ts"]) # Timestep of last snowfall (days) + OUT["ys"] = np.full((gpsum,), 500.0) # Annual snowfall (mm water equivalent) + OUT["subZ"] = np.full((gpsum, nl), grid["max_subZ"]) # Vertical layer depths (m) + OUT["alb_snow"] = np.full((gpsum,), C["alb_fresh"]) # Snow albedo + OUT["snowmass"] = np.zeros((gpsum,)) # Snow mass (m water equivalent) + + if grid.get("doubledepth", False): # Sets layer thicknesses when 'double depth' is active + mask_indices = np.where(grid["mask"] == 1) + split = grid["split"] for n, split_start in enumerate(split[:-1]): - depth_value = (2.0 ** n) * grid['max_subZ'] - OUT['subZ'][mask_indices[0], split_start:split[n + 1]] = depth_value + depth_value = (2.0**n) * grid["max_subZ"] + OUT["subZ"][mask_indices[0], split_start : split[n + 1]] = depth_value - final_depth_value = (2.0 ** len(split)) * grid['max_subZ'] - OUT['subZ'][mask_indices[0], split[-1]:] = final_depth_value + final_depth_value = (2.0 ** len(split)) * grid["max_subZ"] + OUT["subZ"][mask_indices[0], split[-1] :] = final_depth_value ###################################################### # Declare non-initialized variables in `OUT` ###################################################### - OUT['smb_cumulative'] = np.zeros((gpsum,)) # Cumulative climatic mass balance (m w.e.) - OUT['smb'] = np.zeros((gpsum,)) # Climatic mass balance (m w.e.) - OUT['subK'] = np.zeros((gpsum, nl)) # Vertical conductivity (m2 s-1) - OUT['subCeff'] = np.zeros((gpsum, nl)) # Vertical effective heat capacity (J m-3 K) - OUT['subWvol'] = np.zeros((gpsum, nl)) # Vertical volumetric water content (fraction) - OUT['surfH'] = np.zeros((gpsum,)) # Surface height (m) - OUT['Dfreshsnow'] = np.zeros((gpsum,)) # Fresh snow density (kg m-3) - OUT['tstar'] = np.zeros((gpsum,)) # Albedo decay timescale (days) - OUT['runoff_irr_deep_mean'] = np.zeros((gpsum,)) # Runoff irreducible water (m w.e.) - OUT['albedo'] = np.full((gpsum,), C['alb_ice']) # Albedo (fraction) + OUT["smb_cumulative"] = np.zeros((gpsum,)) # Cumulative climatic mass balance (m w.e.) + OUT["smb"] = np.zeros((gpsum,)) # Climatic mass balance (m w.e.) + OUT["subK"] = np.zeros((gpsum, nl)) # Vertical conductivity (m2 s-1) + OUT["subCeff"] = np.zeros((gpsum, nl)) # Vertical effective heat capacity (J m-3 K) + OUT["subWvol"] = np.zeros((gpsum, nl)) # Vertical volumetric water content (fraction) + OUT["surfH"] = np.zeros((gpsum,)) # Surface height (m) + OUT["Dfreshsnow"] = np.zeros((gpsum,)) # Fresh snow density (kg m-3) + OUT["tstar"] = np.zeros((gpsum,)) # Albedo decay timescale (days) + OUT["runoff_irr_deep_mean"] = np.zeros((gpsum,)) # Runoff irreducible water (m w.e.) + OUT["albedo"] = np.full((gpsum,), C["alb_ice"]) # Albedo (fraction) ###################################################### # Initialize `IN` (model input variables) ###################################################### - IN['T'] = np.zeros((gpsum,)) # Air temperature (K) - IN['P'] = np.zeros((gpsum,)) # Precipitation (m w.e.) - IN['snow'] = np.zeros((gpsum,)) # Snow precipitation (m w.e.) - IN['rain'] = np.zeros((gpsum,)) # Rain precipitation (m w.e.) - IN['yearsnow'] = np.zeros((gpsum, nl)) # Annual snow precipitation (m w.e.) - IN['logyearsnow'] = np.zeros((gpsum, nl)) # Log annual snow precipitation (m w.e.) - IN['C'] = np.zeros((gpsum,)) # Cloud cover fraction - IN['WS'] = np.zeros((gpsum,)) # Wind speed (m s-1) - IN['RH'] = np.zeros((gpsum,)) # Relative humidity (fraction) - IN['q'] = np.zeros((gpsum,)) # Specific humidity (g kg-1) - IN['VP'] = np.zeros((gpsum,)) # Vapor pressure (Pa) - IN['Dair'] = np.zeros((gpsum,)) # Air density (kg m-3) - IN['Pres'] = np.zeros((gpsum,)) # Air pressure (Pa) - IN['SWin'] = np.zeros((gpsum,)) # Incoming SW radiation (W m-2) - IN['LWin'] = np.zeros((gpsum,)) # Incoming LW radiation (W m-2) - + IN["T"] = np.zeros((gpsum,)) # Air temperature (K) + IN["P"] = np.zeros((gpsum,)) # Precipitation (m w.e.) + IN["snow"] = np.zeros((gpsum,)) # Snow precipitation (m w.e.) + IN["rain"] = np.zeros((gpsum,)) # Rain precipitation (m w.e.) + IN["yearsnow"] = np.zeros((gpsum, nl)) # Annual snow precipitation (m w.e.) + IN["logyearsnow"] = np.zeros((gpsum, nl)) # Log annual snow precipitation (m w.e.) + IN["C"] = np.zeros((gpsum,)) # Cloud cover fraction + IN["WS"] = np.zeros((gpsum,)) # Wind speed (m s-1) + IN["RH"] = np.zeros((gpsum,)) # Relative humidity (fraction) + IN["q"] = np.zeros((gpsum,)) # Specific humidity (g kg-1) + IN["VP"] = np.zeros((gpsum,)) # Vapor pressure (Pa) + IN["Dair"] = np.zeros((gpsum,)) # Air density (kg m-3) + IN["Pres"] = np.zeros((gpsum,)) # Air pressure (Pa) + IN["SWin"] = np.zeros((gpsum,)) # Incoming SW radiation (W m-2) + IN["LWin"] = np.zeros((gpsum,)) # Incoming LW radiation (W m-2) OUTFILE = {} # Placeholder for output to be saved to files diff --git a/src/ebfm/LOOP_EBM.py b/src/ebfm/LOOP_EBM.py index 4ff8ed11..4937afab 100644 --- a/src/ebfm/LOOP_EBM.py +++ b/src/ebfm/LOOP_EBM.py @@ -1,10 +1,11 @@ import numpy as np -from ebfm import LOOP_EBM_SHF,LOOP_EBM_GHF, LOOP_EBM_LHF, LOOP_EBM_LWin, LOOP_EBM_LWout, LOOP_EBM_SWin +from ebfm import LOOP_EBM_SHF, LOOP_EBM_GHF, LOOP_EBM_LHF, LOOP_EBM_LWin, LOOP_EBM_LWout, LOOP_EBM_SWin from ebfm import LOOP_EBM_SWout, LOOP_EBM_insolation from coupling import Coupler + def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: """ Surface Energy Balance Model: Calculates heat fluxes, surface temperature, melt rates, and moisture fluxes. @@ -30,23 +31,22 @@ def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: # TODO: better do this before calling LOOP_EBM.main if cpl.couple_to_icon_atmo: - LWin = IN['LWin'] + LWin = IN["LWin"] else: LWin = LOOP_EBM_LWin.main(C, IN) SWout, OUT = LOOP_EBM_SWout.main(C, time2, OUT, SWin) GHF_k = 0.138 - 1.01e-3 * OUT["subD"] + 3.233e-6 * OUT["subD"] ** 2 - GHF_C = ( - (GHF_k[:, 0] * OUT["subZ"][:, 0] + 0.5 * GHF_k[:, 1] * OUT["subZ"][:, 1]) - / (OUT["subZ"][:, 0] + 0.5 * OUT["subZ"][:, 1]) ** 2 - ) + GHF_C = (GHF_k[:, 0] * OUT["subZ"][:, 0] + 0.5 * GHF_k[:, 1] * OUT["subZ"][:, 1]) / ( + OUT["subZ"][:, 0] + 0.5 * OUT["subZ"][:, 1] + ) ** 2 # Precompute reusable constant arrays - condition_mask = np.ones(grid['gpsum'], dtype=bool) + condition_mask = np.ones(grid["gpsum"], dtype=bool) # Set initial temperature range - Tlow = OUT['Tsurf'] - 40.0 - Thigh = OUT['Tsurf'] + 40.0 + Tlow = OUT["Tsurf"] - 40.0 + Thigh = OUT["Tsurf"] + 40.0 dT = Thigh - Tlow for c in range(20): @@ -56,20 +56,24 @@ def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: # Surface energy balance at Tlow ebal_Tlow = ( - SWin - SWout + LWin - - LOOP_EBM_LWout.main(C, Tlow) - + LOOP_EBM_LHF.main(C, Tlow, IN, condition_mask) - + LOOP_EBM_SHF.main(C, Tlow, IN, condition_mask) - + LOOP_EBM_GHF.main(Tlow, OUT, condition_mask, GHF_k, GHF_C) + SWin + - SWout + + LWin + - LOOP_EBM_LWout.main(C, Tlow) + + LOOP_EBM_LHF.main(C, Tlow, IN, condition_mask) + + LOOP_EBM_SHF.main(C, Tlow, IN, condition_mask) + + LOOP_EBM_GHF.main(Tlow, OUT, condition_mask, GHF_k, GHF_C) ) # Surface energy balance at Tmid ebal_Tmid = ( - SWin - SWout + LWin - - LOOP_EBM_LWout.main(C, Tmid) - + LOOP_EBM_LHF.main(C, Tmid, IN, condition_mask) - + LOOP_EBM_SHF.main(C, Tmid, IN, condition_mask) - + LOOP_EBM_GHF.main(Tmid, OUT, condition_mask, GHF_k, GHF_C) + SWin + - SWout + + LWin + - LOOP_EBM_LWout.main(C, Tmid) + + LOOP_EBM_LHF.main(C, Tmid, IN, condition_mask) + + LOOP_EBM_SHF.main(C, Tmid, IN, condition_mask) + + LOOP_EBM_GHF.main(Tmid, OUT, condition_mask, GHF_k, GHF_C) ) # Update temperature range based on energy balance sign @@ -78,7 +82,7 @@ def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: Tlow[~cond_EB] = Tmid[~cond_EB] # Stop if temperature change is below the predefined limit - if np.max(dT) < C['dTacc']: + if np.max(dT) < C["dTacc"]: break if c == 19: @@ -89,8 +93,8 @@ def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: ########################################################### # Ensure surface temperature does not exceed the melting point - Tmid[np.abs(Tmid - C['T0']) < C['dTacc']] -= C['dTacc'] - Tmid = np.minimum(Tmid, C['T0']) + Tmid[np.abs(Tmid - C["T0"]) < C["dTacc"]] -= C["dTacc"] + Tmid = np.minimum(Tmid, C["T0"]) LWout = LOOP_EBM_LWout.main(C, Tmid) LHF = LOOP_EBM_LHF.main(C, Tmid, IN, condition_mask) @@ -98,38 +102,38 @@ def main(C, OUT, IN, time2, grid, cpl: Coupler) -> dict: GHF = LOOP_EBM_GHF.main(Tmid, OUT, condition_mask, GHF_k, GHF_C) Emelt = SWin - SWout + LWin - LWout + SHF + LHF + GHF - Emelt[Tmid < C['T0']] = 0.0 + Emelt[Tmid < C["T0"]] = 0.0 - OUT['melt'] = C['dayseconds'] * time2['dt'] * Emelt / C['Lm'] / 1e3 + OUT["melt"] = C["dayseconds"] * time2["dt"] * Emelt / C["Lm"] / 1e3 ########################################################### # MOISTURE FLUXES ########################################################### - OUT['moist_deposition'] = (C['dayseconds'] * time2['dt'] * LHF / C['Ls'] / 1e3 * (Tmid < C['T0']) * (LHF > 0)) - OUT['moist_condensation'] = (C['dayseconds'] * time2['dt'] * LHF / C['Lv'] / 1e3 * (Tmid >= C['T0']) * (LHF > 0)) - OUT['moist_sublimation'] = (-C['dayseconds'] * time2['dt'] * LHF / C['Ls'] / 1e3 * (Tmid < C['T0']) * (LHF < 0)) - OUT['moist_evaporation'] = (-C['dayseconds'] * time2['dt'] * LHF / C['Lv'] / 1e3 * (Tmid >= C['T0']) * (LHF < 0)) + OUT["moist_deposition"] = C["dayseconds"] * time2["dt"] * LHF / C["Ls"] / 1e3 * (Tmid < C["T0"]) * (LHF > 0) + OUT["moist_condensation"] = C["dayseconds"] * time2["dt"] * LHF / C["Lv"] / 1e3 * (Tmid >= C["T0"]) * (LHF > 0) + OUT["moist_sublimation"] = -C["dayseconds"] * time2["dt"] * LHF / C["Ls"] / 1e3 * (Tmid < C["T0"]) * (LHF < 0) + OUT["moist_evaporation"] = -C["dayseconds"] * time2["dt"] * LHF / C["Lv"] / 1e3 * (Tmid >= C["T0"]) * (LHF < 0) ########################################################### # AVOID EVAPORATION OF ABSENT MELT ########################################################### - max_evap = OUT['melt'] - OUT['moist_evaporation'] = np.minimum(OUT['moist_evaporation'], max_evap) + max_evap = OUT["melt"] + OUT["moist_evaporation"] = np.minimum(OUT["moist_evaporation"], max_evap) ########################################################### # STORE RELEVANT VARIABLES IN OUT ########################################################### - OUT['Tsurf'] = Tmid - OUT['Emelt'] = Emelt - OUT['LHF'] = LHF - OUT['SWin'] = SWin - OUT['SWout'] = SWout - OUT['LWin'] = LWin - OUT['LWout'] = LWout - OUT['SHF'] = SHF - OUT['GHF'] = GHF + OUT["Tsurf"] = Tmid + OUT["Emelt"] = Emelt + OUT["LHF"] = LHF + OUT["SWin"] = SWin + OUT["SWout"] = SWout + OUT["LWin"] = LWin + OUT["LWout"] = LWout + OUT["SHF"] = SHF + OUT["GHF"] = GHF return OUT diff --git a/src/ebfm/LOOP_EBM_LHF.py b/src/ebfm/LOOP_EBM_LHF.py index 4eea5adf..651641ce 100644 --- a/src/ebfm/LOOP_EBM_LHF.py +++ b/src/ebfm/LOOP_EBM_LHF.py @@ -36,27 +36,24 @@ def main(C, Tsurf, IN, cond): # Turbulent Exchange Coefficient ########################################################### C_kat = np.maximum( - C["k_turb"] * (IN["T"][cond] - Tsurf) * np.sqrt(C["g"] / (C["T0"] * (IN["Theta_lapse"] * C["Pr"]))), - 0 + C["k_turb"] * (IN["T"][cond] - Tsurf) * np.sqrt(C["g"] / (C["T0"] * (IN["Theta_lapse"] * C["Pr"]))), 0 ) C_turb = 0.5 * (C["turb"] + C_kat) ########################################################### # Saturation Vapor Pressure (Clausius-Clapeyron) ########################################################### - VPsurf = ( - C["VP0"] * np.exp(C["Ls"] / C["Rv"] / 273.15 * (1 - 273.15 / Tsurf)) * (Tsurf < 273.15) - + C["VP0"] * np.exp(C["Lv"] / C["Rv"] / 273.15 * (1 - 273.15 / Tsurf)) * (Tsurf >= 273.15) - ) + VPsurf = C["VP0"] * np.exp(C["Ls"] / C["Rv"] / 273.15 * (1 - 273.15 / Tsurf)) * (Tsurf < 273.15) + C[ + "VP0" + ] * np.exp(C["Lv"] / C["Rv"] / 273.15 * (1 - 273.15 / Tsurf)) * (Tsurf >= 273.15) ########################################################### # Latent Heat Flux (bulk equations) ########################################################### - LHF = ( - C["eps"] * IN["Dair"][cond] * C["Ls"] * C_turb * (IN["VP"][cond] - VPsurf) / IN["Pres"][cond] * ( - Tsurf < 273.15) - + C["eps"] * IN["Dair"][cond] * C["Lv"] * C_turb * (IN["VP"][cond] - VPsurf) / IN["Pres"][cond] * ( - Tsurf >= 273.15) + LHF = C["eps"] * IN["Dair"][cond] * C["Ls"] * C_turb * (IN["VP"][cond] - VPsurf) / IN["Pres"][cond] * ( + Tsurf < 273.15 + ) + C["eps"] * IN["Dair"][cond] * C["Lv"] * C_turb * (IN["VP"][cond] - VPsurf) / IN["Pres"][cond] * ( + Tsurf >= 273.15 ) return LHF diff --git a/src/ebfm/LOOP_EBM_LWin.py b/src/ebfm/LOOP_EBM_LWin.py index be731232..17f82caa 100644 --- a/src/ebfm/LOOP_EBM_LWin.py +++ b/src/ebfm/LOOP_EBM_LWin.py @@ -5,7 +5,7 @@ def main(C, IN): """ Computes incoming longwave radiation (LWin) based on atmospheric emissivity. SOURCE: Konzelmann et al. (1994) - + Parameters: C (dict): Constants including coefficients for emissivity and Stefan-Boltzmann constant. - C["b"]: Coefficient for clear-sky emissivity. diff --git a/src/ebfm/LOOP_EBM_LWout.py b/src/ebfm/LOOP_EBM_LWout.py index cb4e0968..d97b3ad0 100644 --- a/src/ebfm/LOOP_EBM_LWout.py +++ b/src/ebfm/LOOP_EBM_LWout.py @@ -17,5 +17,5 @@ def main(C, Tsurf): ########################################################### # Blackbody Emission of Thermal Radiation ########################################################### - LWout = C["boltz"] * Tsurf ** 4 + LWout = C["boltz"] * Tsurf**4 return LWout diff --git a/src/ebfm/LOOP_EBM_SHF.py b/src/ebfm/LOOP_EBM_SHF.py index e6fab57f..87e9e8c9 100644 --- a/src/ebfm/LOOP_EBM_SHF.py +++ b/src/ebfm/LOOP_EBM_SHF.py @@ -28,8 +28,7 @@ def main(C, Tsurf, IN, cond): # Turbulent Exchange Coefficient ########################################################### C_kat = np.maximum( - C["k_turb"] * (IN["T"][cond] - Tsurf) * np.sqrt(C["g"] / (C["T0"] * (IN["Theta_lapse"] * C["Pr"]))), - 0 + C["k_turb"] * (IN["T"][cond] - Tsurf) * np.sqrt(C["g"] / (C["T0"] * (IN["Theta_lapse"] * C["Pr"]))), 0 ) C_turb = 0.5 * (C["turb"] + C_kat) diff --git a/src/ebfm/LOOP_EBM_SWin.py b/src/ebfm/LOOP_EBM_SWin.py index c3077f86..2a6da964 100644 --- a/src/ebfm/LOOP_EBM_SWin.py +++ b/src/ebfm/LOOP_EBM_SWin.py @@ -20,35 +20,33 @@ def main(C, OUT, IN, grid, cpl: Coupler) -> tuple[np.ndarray, dict]: """ if cpl.couple_to_icon_atmo: - SWin_diff = (0.8 - 0.65 * (1 - IN["C"])) * IN['SWin'] - SWin_dir = (0.2 + 0.65 * (1 - IN["C"])) * (1 - OUT["shade"]) * IN['SWin'] + SWin_diff = (0.8 - 0.65 * (1 - IN["C"])) * IN["SWin"] + SWin_dir = (0.2 + 0.65 * (1 - IN["C"])) * (1 - OUT["shade"]) * IN["SWin"] SWin = SWin_dir + SWin_diff else: ########################################################### # Top of the atmosphere radiation ########################################################### - lat_rad = OUT['lat_rad'] - d_rad = OUT['d_rad'] - h_rad = OUT['h_rad'] + lat_rad = OUT["lat_rad"] + d_rad = OUT["d_rad"] + h_rad = OUT["h_rad"] # On a flat surface OUT["TOAflat"] = OUT["I0"] * ( - np.sin(lat_rad) * np.sin(d_rad) + - np.cos(lat_rad) * np.cos(d_rad) * np.cos(h_rad) + np.sin(lat_rad) * np.sin(d_rad) + np.cos(lat_rad) * np.cos(d_rad) * np.cos(h_rad) ) # SOURCE: Iqbal (1983) # On an inclined surface slopebeta = grid["slope_beta"] # Slope angle slopegamma = grid["slope_gamma"] # Slope azimuth angle OUT["TOA"] = OUT["I0"] * ( - (np.sin(lat_rad) * np.cos(slopebeta) - np.cos(lat_rad) * np.sin(slopebeta) * np.cos( - slopegamma)) * np.sin( - d_rad) + - (np.cos(lat_rad) * np.cos(slopebeta) + np.sin(lat_rad) * np.sin(slopebeta) * np.cos( - slopegamma)) * np.cos( - d_rad) * np.cos(h_rad) + - np.cos(d_rad) * np.sin(slopebeta) * np.sin(slopegamma) * np.sin(h_rad) + (np.sin(lat_rad) * np.cos(slopebeta) - np.cos(lat_rad) * np.sin(slopebeta) * np.cos(slopegamma)) + * np.sin(d_rad) + + (np.cos(lat_rad) * np.cos(slopebeta) + np.sin(lat_rad) * np.sin(slopebeta) * np.cos(slopegamma)) + * np.cos(d_rad) + * np.cos(h_rad) + + np.cos(d_rad) * np.sin(slopebeta) * np.sin(slopegamma) * np.sin(h_rad) ) OUT["TOA"] = np.maximum(OUT["TOA"], 0.0) ########################################################### @@ -60,21 +58,22 @@ def main(C, OUT, IN, grid, cpl: Coupler) -> tuple[np.ndarray, dict]: OUT["TOAdiff"] = (0.8 - 0.65 * (1 - IN["C"])) * OUT["TOA"] OUT["TOAshade"] = OUT["TOAdir"] + OUT["TOAdiff"] - ########################################################### # Atmospheric transmissivity ########################################################### # Transmissivity after gaseous absorption / scattering - m = 35.0 * (IN["Pres"] / C["Pref"]) * ( - 1224.0 * (OUT["TOAflat"] / OUT["I0"]) ** 2 + 1.0) ** -0.5 # SOURCE: Meyers and Dale (1983) + m = ( + 35.0 * (IN["Pres"] / C["Pref"]) * (1224.0 * (OUT["TOAflat"] / OUT["I0"]) ** 2 + 1.0) ** -0.5 + ) # SOURCE: Meyers and Dale (1983) t_rg = 1.021 - 0.084 * np.sqrt( - m * (949.0 * (IN["Pres"] / 1e3) * 1e-5 + 0.051)) # SOURCE: Atwater and Brown Jr (1974) + m * (949.0 * (IN["Pres"] / 1e3) * 1e-5 + 0.051) + ) # SOURCE: Atwater and Brown Jr (1974) # Transmissivity after water vapor absorption - temp_dew_kelvin = (1 / 273.15 - (C["Rv"] / C["Ls"]) * np.log( - IN["q"] * IN["Pres"] / (C["eps"] * C["VP0"]) - )) ** -1 # SOURCE: McDonald (1960) + temp_dew_kelvin = ( + 1 / 273.15 - (C["Rv"] / C["Ls"]) * np.log(IN["q"] * IN["Pres"] / (C["eps"] * C["VP0"])) + ) ** -1 # SOURCE: McDonald (1960) temp_dew_fahr = 32.0 + 1.8 * (temp_dew_kelvin - 273.15) # Assign lambda based on latitude region diff --git a/src/ebfm/LOOP_EBM_SWout.py b/src/ebfm/LOOP_EBM_SWout.py index fabf592a..2edc9c7f 100644 --- a/src/ebfm/LOOP_EBM_SWout.py +++ b/src/ebfm/LOOP_EBM_SWout.py @@ -31,16 +31,14 @@ def main(C, time, OUT, SWin): OUT["tstar"][(OUT["Tsurf"] == C["T0"]) & snow_cond] = C["tstar_wet"] snow_below_T0_cond = (OUT["Tsurf"] < C["T0"]) & snow_cond OUT["tstar"][snow_below_T0_cond] = ( - C["tstar_dry"] + np.minimum(C["T0"] - OUT["Tsurf"][snow_below_T0_cond], 10) * C["tstar_K"] + C["tstar_dry"] + np.minimum(C["T0"] - OUT["Tsurf"][snow_below_T0_cond], 10) * C["tstar_K"] ) alb_snow_decay_cond = snow_cond & (OUT["timelastsnow"] < time["TCUR"]) OUT["alb_snow"][alb_snow_decay_cond] -= ( - np.maximum( - OUT["alb_snow"][alb_snow_decay_cond] - C["alb_firn"], 0.0 - ) - / OUT["tstar"][alb_snow_decay_cond] - * time["dt"] + np.maximum(OUT["alb_snow"][alb_snow_decay_cond] - C["alb_firn"], 0.0) + / OUT["tstar"][alb_snow_decay_cond] + * time["dt"] ) fresh_snow_cond = (OUT["timelastsnow"] == time["TCUR"]) | ice_cond diff --git a/src/ebfm/LOOP_EBM_insolation.py b/src/ebfm/LOOP_EBM_insolation.py index 18cd1fc4..4243ed15 100644 --- a/src/ebfm/LOOP_EBM_insolation.py +++ b/src/ebfm/LOOP_EBM_insolation.py @@ -24,25 +24,29 @@ def main(grid, time2, OUT): # Time as a fraction of a year (in radians) start_of_year = datetime.datetime(time2["TCUR"].year, 1, 1) # January 1 of TCUR's year day_of_year = (time2["TCUR"] - start_of_year).days + ( - time2["TCUR"] - start_of_year).seconds / 86400 # Fractional days + time2["TCUR"] - start_of_year + ).seconds / 86400 # Fractional days trad = 2 * np.pi * day_of_year / 365.242 # Top-of-atmosphere radiation on a surface normal to incident beam OUT["I0"] = 1353.0 * (1.0 + 0.034 * np.cos(trad)) # SOURCE: Meyers and Dale (1983) # Solar declination (d) using an approximation formula - d = (0.322003 - - 22.971 * np.cos(trad) - - 0.357898 * np.cos(2.0 * trad) - - 0.14398 * np.cos(3.0 * trad) - + 3.94638 * np.sin(trad) - + 0.019334 * np.sin(2.0 * trad) - + 0.05928 * np.sin(3.0 * trad)) # SOURCE: Iqbal (1983) + d = ( + 0.322003 + - 22.971 * np.cos(trad) + - 0.357898 * np.cos(2.0 * trad) + - 0.14398 * np.cos(3.0 * trad) + + 3.94638 * np.sin(trad) + + 0.019334 * np.sin(2.0 * trad) + + 0.05928 * np.sin(3.0 * trad) + ) # SOURCE: Iqbal (1983) # Solar hour angle (h) B = 360 / 365 * (day_of_year - 81) - Tcor_ecc = 9.87 * np.sin(2.0 * np.radians(B)) - 7.53 * np.cos(np.radians(B)) - 1.5 * np.sin( - np.radians(B)) # Correction for eccentricity + Tcor_ecc = ( + 9.87 * np.sin(2.0 * np.radians(B)) - 7.53 * np.cos(np.radians(B)) - 1.5 * np.sin(np.radians(B)) + ) # Correction for eccentricity Tcor_lon = 4 * (grid["lon"] - 15 * time2["dT_UTC"]) # Correction for longitude within time-zone Tcor = Tcor_ecc + Tcor_lon LST = time2["TCUR"].hour + time2["TCUR"].minute / 60 + Tcor / 60 # Local Solar Time @@ -52,21 +56,20 @@ def main(grid, time2, OUT): d_rad = np.radians(d) h_rad = np.radians(h) - OUT['lat_rad'] = lat_rad - OUT['d_rad'] = d_rad - OUT['h_rad'] = h_rad + OUT["lat_rad"] = lat_rad + OUT["d_rad"] = d_rad + OUT["h_rad"] = h_rad ########################################################### # Shading by the surrounding topography ########################################################### elevationangle = np.arcsin( - np.sin(lat_rad) * np.sin(d_rad) + - np.cos(lat_rad) * np.cos(d_rad) * np.cos(h_rad) + np.sin(lat_rad) * np.sin(d_rad) + np.cos(lat_rad) * np.cos(d_rad) * np.cos(h_rad) ) # SOURCE: Iqbal (1983) - if grid['input_type'] is GridInputType.MATLAB: # only use shading routine in uncoupled runs on a regular grid - yl = grid['Ly'] - xl = grid['Lx'] + if grid["input_type"] is GridInputType.MATLAB: # only use shading routine in uncoupled runs on a regular grid + yl = grid["Ly"] + xl = grid["Lx"] # Azimuth calculation cos_elevation = np.cos(elevationangle) @@ -77,7 +80,7 @@ def main(grid, time2, OUT): ), -np.arccos( (np.cos(h_rad) * np.cos(d_rad) * np.sin(lat_rad) - np.sin(d_rad) * np.cos(lat_rad)) / cos_elevation - ) + ), ) # Directional gradients based on azimuth @@ -87,39 +90,42 @@ def main(grid, time2, OUT): ddx[azimuth <= -0.75 * np.pi] = -np.tan(np.pi + azimuth[azimuth <= -0.75 * np.pi]) ddx[(azimuth > -0.75 * np.pi) & (azimuth <= -0.25 * np.pi)] = -1 ddx[(azimuth > -0.25 * np.pi) & (azimuth <= 0.25 * np.pi)] = np.tan( - azimuth[(azimuth > -0.25 * np.pi) & (azimuth <= 0.25 * np.pi)]) + azimuth[(azimuth > -0.25 * np.pi) & (azimuth <= 0.25 * np.pi)] + ) ddx[(azimuth > 0.25 * np.pi) & (azimuth <= 0.75 * np.pi)] = 1 ddx[azimuth > 0.75 * np.pi] = np.tan(np.pi - azimuth[azimuth > 0.75 * np.pi]) ddy[azimuth <= -0.75 * np.pi] = 1 ddy[(azimuth > -0.75 * np.pi) & (azimuth <= -0.25 * np.pi)] = -np.tan( - np.pi / 2 + azimuth[(azimuth > -0.75 * np.pi) & (azimuth <= -0.25 * np.pi)]) + np.pi / 2 + azimuth[(azimuth > -0.75 * np.pi) & (azimuth <= -0.25 * np.pi)] + ) ddy[(azimuth > -0.25 * np.pi) & (azimuth <= 0.25 * np.pi)] = -1 ddy[(azimuth > 0.25 * np.pi) & (azimuth <= 0.75 * np.pi)] = -np.tan( - np.pi / 2 - azimuth[(azimuth > 0.25 * np.pi) & (azimuth <= 0.75 * np.pi)]) + np.pi / 2 - azimuth[(azimuth > 0.25 * np.pi) & (azimuth <= 0.75 * np.pi)] + ) ddy[azimuth > 0.75 * np.pi] = 1 # Initialize the shade_2D array mask = grid["mask_2D"] # Shading routine - z_flat = grid['z_2D'].ravel() - z_current = z_flat[grid['ind']] - shade = np.zeros(grid['gpsum'], dtype=int) + z_flat = grid["z_2D"].ravel() + z_current = z_flat[grid["ind"]] + shade = np.zeros(grid["gpsum"], dtype=int) max_count = 200 # Precompute grid indices and deltas ahead of time delta_y, delta_x = np.round(ddx).astype(int), np.round(ddy).astype(int) # Precompute flattened grid distances - ddx_scaled = ddx * grid['dx'] - ddy_scaled = ddy * grid['dx'] - distance_grid = np.sqrt(ddx_scaled ** 2 + ddy_scaled ** 2) + ddx_scaled = ddx * grid["dx"] + ddy_scaled = ddy * grid["dx"] + distance_grid = np.sqrt(ddx_scaled**2 + ddy_scaled**2) distance_scaled = np.outer(np.arange(1, max_count + 1), distance_grid) for count in range(1, max_count + 1): - kk = np.clip(grid['yind'] + delta_y * count, 0, yl - 1) - ll = np.clip(grid['xind'] + delta_x * count, 0, xl - 1) + kk = np.clip(grid["yind"] + delta_y * count, 0, yl - 1) + ll = np.clip(grid["xind"] + delta_x * count, 0, xl - 1) ind_kkll = ll * yl + kk elev_diff = z_flat[ind_kkll] - z_current @@ -134,12 +140,12 @@ def main(grid, time2, OUT): # Build the final shade 2D array shade_2D = np.ones((xl, yl), dtype=int) - shade_2D.flat[grid['ind']] = shade + shade_2D.flat[grid["ind"]] = shade OUT["shade"] = shade_2D.flatten()[grid["mask_2D"].flatten() == 1] else: # TODO: Shading based on lookup tables from Elmer/Ice # For now: assume flat surface and no shading by surrounding terrain - OUT["shade"] = np.zeros_like(grid['x']) - OUT["shade"][elevationangle<0] = 1 + OUT["shade"] = np.zeros_like(grid["x"]) + OUT["shade"][elevationangle < 0] = 1 return OUT diff --git a/src/ebfm/LOOP_SNOW.py b/src/ebfm/LOOP_SNOW.py index ac08a8c9..b95b0efc 100644 --- a/src/ebfm/LOOP_SNOW.py +++ b/src/ebfm/LOOP_SNOW.py @@ -21,35 +21,36 @@ def snowfall_and_deposition(): Calculate snowfall and deposition and shift vertical grid accordingly """ - max_subZ = grid['max_subZ'] - gpsum = grid['gpsum'] - nl = grid['nl'] + max_subZ = grid["max_subZ"] + gpsum = grid["gpsum"] + nl = grid["nl"] # Fresh snow density calculations - if phys['snow_compaction'] == 'firn+snow': - OUT['Dfreshsnow_T'] = np.zeros_like(IN['T']) + if phys["snow_compaction"] == "firn+snow": + OUT["Dfreshsnow_T"] = np.zeros_like(IN["T"]) - temp_above = IN['T'] > C['T0'] + 2 - temp_between = (IN['T'] <= C['T0'] + 2) & (IN['T'] > C['T0'] - 15) - temp_below = IN['T'] <= C['T0'] - 15 + temp_above = IN["T"] > C["T0"] + 2 + temp_between = (IN["T"] <= C["T0"] + 2) & (IN["T"] > C["T0"] - 15) + temp_below = IN["T"] <= C["T0"] - 15 - OUT['Dfreshsnow_T'][temp_above] = 50 + 1.7 * 17 ** (3 / 2) - OUT['Dfreshsnow_T'][temp_between] = 50 + 1.7 * (IN['T'][temp_between] - C['T0'] + 15) ** (3 / 2) - OUT['Dfreshsnow_T'][temp_below] = -3.8328 * (IN['T'][temp_below] - C['T0']) - \ - 0.0333 * (IN['T'][temp_below] - C['T0']) ** 2 + OUT["Dfreshsnow_T"][temp_above] = 50 + 1.7 * 17 ** (3 / 2) + OUT["Dfreshsnow_T"][temp_between] = 50 + 1.7 * (IN["T"][temp_between] - C["T0"] + 15) ** (3 / 2) + OUT["Dfreshsnow_T"][temp_below] = ( + -3.8328 * (IN["T"][temp_below] - C["T0"]) - 0.0333 * (IN["T"][temp_below] - C["T0"]) ** 2 + ) - OUT['Dfreshsnow_W'] = 266.86 * (0.5 * (1 + np.tanh(IN['WS'] / 5))) ** 8.8 - OUT['Dfreshsnow'] = OUT['Dfreshsnow_T'] + OUT['Dfreshsnow_W'] + OUT["Dfreshsnow_W"] = 266.86 * (0.5 * (1 + np.tanh(IN["WS"] / 5))) ** 8.8 + OUT["Dfreshsnow"] = OUT["Dfreshsnow_T"] + OUT["Dfreshsnow_W"] else: - OUT['Dfreshsnow'] = np.full_like(IN['snow'], C['Dfreshsnow']) + OUT["Dfreshsnow"] = np.full_like(IN["snow"], C["Dfreshsnow"]) # Update layer depths and properties - shift_snowfall = IN['snow'] * C['Dwater'] / OUT['Dfreshsnow'] - shift_riming = OUT['moist_deposition'] * C['Dwater'] / OUT['Dfreshsnow'] + shift_snowfall = IN["snow"] * C["Dwater"] / OUT["Dfreshsnow"] + shift_riming = OUT["moist_deposition"] * C["Dwater"] / OUT["Dfreshsnow"] shift_tot = shift_snowfall + shift_riming - OUT['surfH'] += shift_tot + OUT["surfH"] += shift_tot - OUT['runoff_irr_deep'] = np.zeros(gpsum) + OUT["runoff_irr_deep"] = np.zeros(gpsum) # Main processing loop: Run until all shifts are handled while np.any(shift_tot > 0): @@ -57,52 +58,54 @@ def snowfall_and_deposition(): shift_tot -= shift # Use references instead of unnecessary .copy() - subT_old = OUT['subT'].copy() - subD_old = OUT['subD'].copy() - subW_old = OUT['subW'].copy() - subZ_old = OUT['subZ'].copy() + subT_old = OUT["subT"].copy() + subD_old = OUT["subD"].copy() + subW_old = OUT["subW"].copy() + subZ_old = OUT["subZ"].copy() # Precompute conditions for better performance is_noshift = subZ_old[:, 0] + shift <= max_subZ is_shift = ~is_noshift # Handle no-shift updates (vectorized) - OUT['subZ'][is_noshift, 0] += shift[is_noshift] - OUT['subT'][is_noshift, 0] = ( - subT_old[is_noshift, 0] * subZ_old[is_noshift, 0] / OUT['subZ'][is_noshift, 0] + - OUT['Tsurf'][is_noshift] * shift[is_noshift] / OUT['subZ'][is_noshift, 0] + OUT["subZ"][is_noshift, 0] += shift[is_noshift] + OUT["subT"][is_noshift, 0] = ( + subT_old[is_noshift, 0] * subZ_old[is_noshift, 0] / OUT["subZ"][is_noshift, 0] + + OUT["Tsurf"][is_noshift] * shift[is_noshift] / OUT["subZ"][is_noshift, 0] ) - OUT['subD'][is_noshift, 0] = ( - subD_old[is_noshift, 0] * subZ_old[is_noshift, 0] / OUT['subZ'][is_noshift, 0] + - OUT['Dfreshsnow'][is_noshift] * shift[is_noshift] / OUT['subZ'][is_noshift, 0] + OUT["subD"][is_noshift, 0] = ( + subD_old[is_noshift, 0] * subZ_old[is_noshift, 0] / OUT["subZ"][is_noshift, 0] + + OUT["Dfreshsnow"][is_noshift] * shift[is_noshift] / OUT["subZ"][is_noshift, 0] ) # Handle shifting updates (vectorized) if np.any(is_shift): - OUT['subZ'][is_shift, 2:nl - 1] = subZ_old[is_shift, 1:nl - 2] - OUT['subT'][is_shift, 2:nl - 1] = subT_old[is_shift, 1:nl - 2] - OUT['subD'][is_shift, 2:nl - 1] = subD_old[is_shift, 1:nl - 2] - OUT['subW'][is_shift, 2:nl - 1] = subW_old[is_shift, 1:nl - 2] - - OUT['subZ'][is_shift, 1] = max_subZ - OUT['subZ'][is_shift, 0] = (subZ_old[is_shift, 0] + shift[is_shift]) - max_subZ - OUT['subT'][is_shift, 1] = ( - subT_old[is_shift, 0] * subZ_old[is_shift, 0] / OUT['subZ'][is_shift, 1] + - OUT['Tsurf'][is_shift] * (OUT['subZ'][is_shift, 1] - subZ_old[is_shift, 0]) / - OUT['subZ'][is_shift, 1] + OUT["subZ"][is_shift, 2 : nl - 1] = subZ_old[is_shift, 1 : nl - 2] + OUT["subT"][is_shift, 2 : nl - 1] = subT_old[is_shift, 1 : nl - 2] + OUT["subD"][is_shift, 2 : nl - 1] = subD_old[is_shift, 1 : nl - 2] + OUT["subW"][is_shift, 2 : nl - 1] = subW_old[is_shift, 1 : nl - 2] + + OUT["subZ"][is_shift, 1] = max_subZ + OUT["subZ"][is_shift, 0] = (subZ_old[is_shift, 0] + shift[is_shift]) - max_subZ + OUT["subT"][is_shift, 1] = ( + subT_old[is_shift, 0] * subZ_old[is_shift, 0] / OUT["subZ"][is_shift, 1] + + OUT["Tsurf"][is_shift] + * (OUT["subZ"][is_shift, 1] - subZ_old[is_shift, 0]) + / OUT["subZ"][is_shift, 1] ) - OUT['subT'][is_shift, 0] = OUT['Tsurf'][is_shift] - OUT['subD'][is_shift, 1] = ( - subD_old[is_shift, 0] * subZ_old[is_shift, 0] / OUT['subZ'][is_shift, 1] + - OUT['Dfreshsnow'][is_shift] * (OUT['subZ'][is_shift, 1] - subZ_old[is_shift, 0]) / - OUT['subZ'][is_shift, 1] + OUT["subT"][is_shift, 0] = OUT["Tsurf"][is_shift] + OUT["subD"][is_shift, 1] = ( + subD_old[is_shift, 0] * subZ_old[is_shift, 0] / OUT["subZ"][is_shift, 1] + + OUT["Dfreshsnow"][is_shift] + * (OUT["subZ"][is_shift, 1] - subZ_old[is_shift, 0]) + / OUT["subZ"][is_shift, 1] ) - OUT['subD'][is_shift, 0] = OUT['Dfreshsnow'][is_shift] - OUT['subW'][is_shift, 1] = subW_old[is_shift, 0] - OUT['subW'][is_shift, 0] = 0.0 + OUT["subD"][is_shift, 0] = OUT["Dfreshsnow"][is_shift] + OUT["subW"][is_shift, 1] = subW_old[is_shift, 0] + OUT["subW"][is_shift, 0] = 0.0 # Update runoff for shifted layers - OUT['runoff_irr_deep'][is_shift] += subW_old[is_shift, nl - 1] + OUT["runoff_irr_deep"][is_shift] += subW_old[is_shift, nl - 1] return True @@ -112,9 +115,9 @@ def melt_sublimation(): """ # Initialize variables - OUT['sumWinit'] = np.sum(OUT['subW'], axis=1) - mass_removed = (OUT['melt'] + OUT['moist_sublimation']) * 1e3 - mass_layer = OUT['subD'] * OUT['subZ'] + OUT["sumWinit"] = np.sum(OUT["subW"], axis=1) + mass_removed = (OUT["melt"] + OUT["moist_sublimation"]) * 1e3 + mass_layer = OUT["subD"] * OUT["subZ"] shift_tot = np.zeros_like(mass_removed) n = 0 @@ -126,59 +129,59 @@ def melt_sublimation(): cond2 = (~cond1) & (mass_removed > 0) # Update for layers fully removed - mass_removed[cond1] -= OUT['subD'][cond1, n - 1] * OUT['subZ'][cond1, n - 1] - shift_tot[cond1] -= OUT['subZ'][cond1, n - 1] + mass_removed[cond1] -= OUT["subD"][cond1, n - 1] * OUT["subZ"][cond1, n - 1] + shift_tot[cond1] -= OUT["subZ"][cond1, n - 1] # Update for layers partially removed - shift_tot[cond2] -= (mass_removed[cond2] / mass_layer[cond2, n - 1]) * OUT['subZ'][cond2, n - 1] + shift_tot[cond2] -= (mass_removed[cond2] / mass_layer[cond2, n - 1]) * OUT["subZ"][cond2, n - 1] mass_removed[cond2] = 0.0 # While there are shifts required while np.any(shift_tot < 0): - shift = np.maximum(shift_tot, -OUT['subZ'][:, 1]) + shift = np.maximum(shift_tot, -OUT["subZ"][:, 1]) shift_tot -= shift - OUT['surfH'] += shift + OUT["surfH"] += shift # Save old values for updates - subT_old = OUT['subT'].copy() - subD_old = OUT['subD'].copy() - subW_old = OUT['subW'].copy() - subZ_old = OUT['subZ'].copy() + subT_old = OUT["subT"].copy() + subD_old = OUT["subD"].copy() + subW_old = OUT["subW"].copy() + subZ_old = OUT["subZ"].copy() # Find no-shift and shift indices i_noshift = np.where(subZ_old[:, 0] + shift > 1e-17) i_shift = np.where(subZ_old[:, 0] + shift <= 1e-17) # Handle the no-shift case - OUT['subZ'][i_noshift, 0] = subZ_old[i_noshift, 0] + shift[i_noshift] - OUT['subT'][i_noshift, 0] = subT_old[i_noshift, 0] - OUT['subD'][i_noshift, 0] = subD_old[i_noshift, 0] - temp = OUT['subZ'][i_noshift, 0] / subZ_old[i_noshift, 0] - OUT['subW'][i_noshift, 0] = subW_old[i_noshift, 0] * temp + OUT["subZ"][i_noshift, 0] = subZ_old[i_noshift, 0] + shift[i_noshift] + OUT["subT"][i_noshift, 0] = subT_old[i_noshift, 0] + OUT["subD"][i_noshift, 0] = subD_old[i_noshift, 0] + temp = OUT["subZ"][i_noshift, 0] / subZ_old[i_noshift, 0] + OUT["subW"][i_noshift, 0] = subW_old[i_noshift, 0] * temp # Handle the shift case - nl = grid['nl'] - OUT['subZ'][i_shift, 1:nl - 2] = subZ_old[i_shift, 2:nl - 1] - OUT['subT'][i_shift, 1:nl - 2] = subT_old[i_shift, 2:nl - 1] - OUT['subD'][i_shift, 1:nl - 2] = subD_old[i_shift, 2:nl - 1] - OUT['subW'][i_shift, 1:nl - 2] = subW_old[i_shift, 2:nl - 1] - - OUT['subZ'][i_shift, 0] = subZ_old[i_shift, 0] + subZ_old[i_shift, 1] + shift[i_shift] - OUT['subT'][i_shift, 0] = subT_old[i_shift, 1] - OUT['subD'][i_shift, 0] = subD_old[i_shift, 1] - temp = OUT['subZ'][i_shift, 0] / subZ_old[i_shift, 1] - OUT['subW'][i_shift, 0] = subW_old[i_shift, 1] * temp - OUT['subT'][i_shift, nl - 1] = subT_old[i_shift, nl - 1] - OUT['subW'][i_shift, nl - 1] = 0.0 + nl = grid["nl"] + OUT["subZ"][i_shift, 1 : nl - 2] = subZ_old[i_shift, 2 : nl - 1] + OUT["subT"][i_shift, 1 : nl - 2] = subT_old[i_shift, 2 : nl - 1] + OUT["subD"][i_shift, 1 : nl - 2] = subD_old[i_shift, 2 : nl - 1] + OUT["subW"][i_shift, 1 : nl - 2] = subW_old[i_shift, 2 : nl - 1] + + OUT["subZ"][i_shift, 0] = subZ_old[i_shift, 0] + subZ_old[i_shift, 1] + shift[i_shift] + OUT["subT"][i_shift, 0] = subT_old[i_shift, 1] + OUT["subD"][i_shift, 0] = subD_old[i_shift, 1] + temp = OUT["subZ"][i_shift, 0] / subZ_old[i_shift, 1] + OUT["subW"][i_shift, 0] = subW_old[i_shift, 1] * temp + OUT["subT"][i_shift, nl - 1] = subT_old[i_shift, nl - 1] + OUT["subW"][i_shift, nl - 1] = 0.0 # Update the deepest layer properties for idx in i_shift: - if grid['doubledepth'] == 1: - OUT['subZ'][idx, nl - 1] = 2.0 ** len(grid['split']) * grid['max_subZ'] + if grid["doubledepth"] == 1: + OUT["subZ"][idx, nl - 1] = 2.0 ** len(grid["split"]) * grid["max_subZ"] else: - OUT['subZ'][idx, nl - 1] = grid['max_subZ'] - OUT['subD'][idx, nl - 1] = subD_old[idx, nl - 1] + OUT["subZ"][idx, nl - 1] = grid["max_subZ"] + OUT["subD"][idx, nl - 1] = subD_old[idx, nl - 1] return True @@ -187,137 +190,140 @@ def compaction(): Calculate snow and firn compaction and update density and layer thickness """ - gpsum, nl = grid['gpsum'], grid['nl'] - Dice, Dfirn = C['Dice'], C['Dfirn'] + gpsum, nl = grid["gpsum"], grid["nl"] + Dice, Dfirn = C["Dice"], C["Dfirn"] - subD_old = OUT['subD'].copy() - subZ_old = OUT['subZ'].copy() + subD_old = OUT["subD"].copy() + subZ_old = OUT["subZ"].copy() mliqmax = np.zeros((gpsum, nl)) - dt_yearfrac = dt / C['yeardays'] - dt_seconds = dt * C['dayseconds'] + dt_yearfrac = dt / C["yeardays"] + dt_seconds = dt * C["dayseconds"] # ------ FIRN COMPACTION ------ # - if phys['snow_compaction'] in ['firn_only', 'firn+snow']: + if phys["snow_compaction"] in ["firn_only", "firn+snow"]: # Pre-compute the logical condition based on the snow compaction type - if phys['snow_compaction'] == 'firn_only': - cond_firn = np.ones_like(OUT['subD'], dtype=bool) # All values are True in 2D + if phys["snow_compaction"] == "firn_only": + cond_firn = np.ones_like(OUT["subD"], dtype=bool) # All values are True in 2D else: # 'firn+snow' - cond_firn = OUT['subD'] >= Dfirn # Results in the same 2D shape as OUT['subD'] + cond_firn = OUT["subD"] >= Dfirn # Results in the same 2D shape as OUT['subD'] # Update annual running average subsurface temperature - OUT['subTmean'] *= (1 - dt_yearfrac) - OUT['subTmean'] += dt_yearfrac * OUT['subT'] + OUT["subTmean"] *= 1 - dt_yearfrac + OUT["subTmean"] += dt_yearfrac * OUT["subT"] # Set gravitational constants - subD_cond = np.where(cond_firn, OUT['subD'], 0) # Using a masked version of subD - logyearsnow_cond = np.where(cond_firn, IN['logyearsnow'], 0) - grav_const = np.zeros_like(OUT['subD']) # Allocation happens here. + subD_cond = np.where(cond_firn, OUT["subD"], 0) # Using a masked version of subD + logyearsnow_cond = np.where(cond_firn, IN["logyearsnow"], 0) + grav_const = np.zeros_like(OUT["subD"]) # Allocation happens here. low_density_mask = cond_firn & (subD_cond < 550) high_density_mask = cond_firn & (subD_cond >= 550) grav_const[low_density_mask] = 0.07 * np.maximum(1.435 - 0.151 * logyearsnow_cond[low_density_mask], 0.25) grav_const[high_density_mask] = 0.03 * np.maximum(2.366 - 0.293 * logyearsnow_cond[high_density_mask], 0.25) # Update firn densities - temp_factor = np.exp(-C['Ec'] / (C['rd'] * OUT['subT']) + C['Eg'] / (C['rd'] * OUT['subTmean'])) - firn_increment = ( - dt_yearfrac * grav_const * IN['yearsnow'] * C['g'] * - (Dice - OUT['subD']) * temp_factor - ) - OUT['subD'][cond_firn] += firn_increment[cond_firn] + temp_factor = np.exp(-C["Ec"] / (C["rd"] * OUT["subT"]) + C["Eg"] / (C["rd"] * OUT["subTmean"])) + firn_increment = dt_yearfrac * grav_const * IN["yearsnow"] * C["g"] * (Dice - OUT["subD"]) * temp_factor + OUT["subD"][cond_firn] += firn_increment[cond_firn] else: raise ValueError("phys.snow_compaction not set correctly!") # ------ SEASONAL SNOW COMPACTION ------ # - if phys['snow_compaction'] == 'firn+snow': + if phys["snow_compaction"] == "firn+snow": # ------ DENSIFICATION BY DESTRUCTIVE METAMORPHISM ------ # # Precompute condition for snow compaction - cond_snow = OUT['subD'] < Dfirn + cond_snow = OUT["subD"] < Dfirn # Constants for snow compaction CC3, CC4 = 2.777e-6, 0.04 # Precompute CC1, CC2, and temp_exp - CC1 = np.exp(-0.046 * np.clip(OUT['subD'] - 175, 0, None)) - CC2 = 1 + (OUT['subW'] != 0) - temp_exp = np.exp(CC4 * (OUT['subT'] - C['T0'])) # + CC1 = np.exp(-0.046 * np.clip(OUT["subD"] - 175, 0, None)) + CC2 = 1 + (OUT["subW"] != 0) + temp_exp = np.exp(CC4 * (OUT["subT"] - C["T0"])) # # Snow densification increment - snow_increment = CC1 * CC2 * CC3 * temp_exp * dt_seconds * OUT['subD'] + snow_increment = CC1 * CC2 * CC3 * temp_exp * dt_seconds * OUT["subD"] # Apply snow increment only to relevant layers - OUT['subD'][cond_snow] += snow_increment[cond_snow] - OUT['subD'][cond_snow] = np.minimum(OUT['subD'][cond_snow], Dice) + OUT["subD"][cond_snow] += snow_increment[cond_snow] + OUT["subD"][cond_snow] = np.minimum(OUT["subD"][cond_snow], Dice) # Store densification by destructive metamorphism - OUT['Dens_destr_metam'] = np.zeros_like(OUT['subD']) - OUT['Dens_destr_metam'][cond_snow] = snow_increment[cond_snow] + OUT["Dens_destr_metam"] = np.zeros_like(OUT["subD"]) + OUT["Dens_destr_metam"][cond_snow] = snow_increment[cond_snow] # ------ DENSIFICATION BY OVERBURDEN PRESSURE ------ # CC5, CC6 = 0.1, 0.023 - CC7 = 4.0 * 7.62237e6 / 250.0 * OUT['subD'] * 1 / (1 + 60 * OUT['subW'] * 1 / (C['Dwater'] * OUT['subZ'])) + CC7 = 4.0 * 7.62237e6 / 250.0 * OUT["subD"] * 1 / (1 + 60 * OUT["subW"] * 1 / (C["Dwater"] * OUT["subZ"])) # Compute load pressure (Psload) - OUT_subD_Z = OUT['subD'] * OUT['subZ'] * C['g'] + OUT_subD_Z = OUT["subD"] * OUT["subZ"] * C["g"] Psload = np.cumsum(OUT_subD_Z, axis=1) - 0.5 * OUT_subD_Z Psload[~cond_snow] = 0 # Compute viscosity (Visc) - temperature_diff = C['T0'] - OUT['subT'] - Visc = CC7 * np.exp(CC5 * temperature_diff + CC6 * OUT['subD']) + temperature_diff = C["T0"] - OUT["subT"] + Visc = CC7 * np.exp(CC5 * temperature_diff + CC6 * OUT["subD"]) Visc[~cond_snow] = 0 # Update densities - OUT['subD'][cond_snow] += dt * C['dayseconds'] * OUT['subD'][cond_snow] * Psload[cond_snow] / Visc[ - cond_snow] - OUT['subD'][cond_snow] = np.minimum(OUT['subD'][cond_snow], C['Dice']) + OUT["subD"][cond_snow] += ( + dt * C["dayseconds"] * OUT["subD"][cond_snow] * Psload[cond_snow] / Visc[cond_snow] + ) + OUT["subD"][cond_snow] = np.minimum(OUT["subD"][cond_snow], C["Dice"]) # Store densification by overburden pressure - OUT['Dens_overb_pres'] = np.zeros_like(OUT['subD']) - OUT['Dens_overb_pres'][cond_snow] = ( - dt * C['dayseconds'] * OUT['subD'][cond_snow] * Psload[cond_snow] / Visc[cond_snow] + OUT["Dens_overb_pres"] = np.zeros_like(OUT["subD"]) + OUT["Dens_overb_pres"][cond_snow] = ( + dt * C["dayseconds"] * OUT["subD"][cond_snow] * Psload[cond_snow] / Visc[cond_snow] ) # ------ DRIFTING SNOW DENSIFICATION ------ # - MO = -0.069 + 0.66 * (1.25 - 0.0042 * (np.maximum(OUT['subD'], 50) - 50)) - SI = -2.868 * np.exp(-0.085 * np.tile(IN['WS'], (nl, 1)).T) + 1 + MO + MO = -0.069 + 0.66 * (1.25 - 0.0042 * (np.maximum(OUT["subD"], 50) - 50)) + SI = -2.868 * np.exp(-0.085 * np.tile(IN["WS"], (nl, 1)).T) + 1 + MO cond_drift = SI > 0 - z_i = np.zeros_like(OUT['subZ']) + z_i = np.zeros_like(OUT["subZ"]) if nl > 1: - z_i[:, 1:] = np.cumsum(OUT['subZ'][:, :-1] * (3.25 - SI[:, :-1]), axis=1) + z_i[:, 1:] = np.cumsum(OUT["subZ"][:, :-1] * (3.25 - SI[:, :-1]), axis=1) gamma_drift = np.maximum(0, SI * np.exp(-z_i / 0.1)) tau = 48 * 2 * 3600 - np.seterr(divide='ignore') + np.seterr(divide="ignore") tau_i = tau / gamma_drift # Update densities - drift_increment = dt_seconds * np.maximum(350 - OUT['subD'], 0) / tau_i - cond_drift_total = cond_drift & (OUT['subD'] < Dfirn) - OUT['subD'][cond_drift_total] += drift_increment[cond_drift_total] - OUT['subD'][cond_drift_total] = np.minimum(OUT['subD'][cond_drift_total], Dice) + drift_increment = dt_seconds * np.maximum(350 - OUT["subD"], 0) / tau_i + cond_drift_total = cond_drift & (OUT["subD"] < Dfirn) + OUT["subD"][cond_drift_total] += drift_increment[cond_drift_total] + OUT["subD"][cond_drift_total] = np.minimum(OUT["subD"][cond_drift_total], Dice) # Store densification by wind shearing - OUT['Dens_drift'] = np.zeros_like(OUT['subD']) - OUT['Dens_drift'][cond_drift_total] = drift_increment[cond_drift_total] + OUT["Dens_drift"] = np.zeros_like(OUT["subD"]) + OUT["Dens_drift"][cond_drift_total] = drift_increment[cond_drift_total] # ------ UPDATE LAYER THICKNESS & SURFACE HEIGHT AFTER COMPACTION ------ # - cond_layers = OUT['subD'] < Dice + cond_layers = OUT["subD"] < Dice # Update layer thickness - OUT['subZ'][cond_layers] = subZ_old[cond_layers] * subD_old[cond_layers] / OUT['subD'][cond_layers] + OUT["subZ"][cond_layers] = subZ_old[cond_layers] * subD_old[cond_layers] / OUT["subD"][cond_layers] # Update irreducible water storage - exp_factor = 0.0143 * np.exp(3.3 * (Dice - OUT['subD'][cond_layers]) / Dice) - mliqmax[cond_layers] = (OUT['subD'][cond_layers] * OUT['subZ'][cond_layers] * exp_factor / - (1 - exp_factor) * 0.05 * - np.minimum(Dice - OUT['subD'][cond_layers], 20)) - OUT['subW'] = np.minimum(mliqmax, OUT['subW']) + exp_factor = 0.0143 * np.exp(3.3 * (Dice - OUT["subD"][cond_layers]) / Dice) + mliqmax[cond_layers] = ( + OUT["subD"][cond_layers] + * OUT["subZ"][cond_layers] + * exp_factor + / (1 - exp_factor) + * 0.05 + * np.minimum(Dice - OUT["subD"][cond_layers], 20) + ) + OUT["subW"] = np.minimum(mliqmax, OUT["subW"]) # Update surface height and runoff - shift = np.sum(OUT['subZ'], axis=1) - np.sum(subZ_old, axis=1) - OUT['surfH'] += shift - OUT['runoff_irr'] = OUT['sumWinit'] - np.sum(OUT['subW'], axis=1) + shift = np.sum(OUT["subZ"], axis=1) - np.sum(subZ_old, axis=1) + OUT["surfH"] += shift + OUT["runoff_irr"] = OUT["sumWinit"] - np.sum(OUT["subW"], axis=1) return True @@ -325,30 +331,27 @@ def heat_conduction(): """ Calculate heat diffusion and update temperatures """ - nl = grid['nl'] - dz1 = (OUT['subZ'][:, 0] + 0.5 * OUT['subZ'][:, 1]) ** 2 - dz2 = 0.5 * (OUT['subZ'][:, 2:] + OUT['subZ'][:, 1:-1]) ** 2 - kk = 0.138 - 1.01e-3 * OUT['subD'] + 3.233e-6 * OUT['subD'] ** 2 # Effective conductivity - c_eff = OUT['subD'] * (152.2 + 7.122 * OUT['subT']) # Effective heat capacity + nl = grid["nl"] + dz1 = (OUT["subZ"][:, 0] + 0.5 * OUT["subZ"][:, 1]) ** 2 + dz2 = 0.5 * (OUT["subZ"][:, 2:] + OUT["subZ"][:, 1:-1]) ** 2 + kk = 0.138 - 1.01e-3 * OUT["subD"] + 3.233e-6 * OUT["subD"] ** 2 # Effective conductivity + c_eff = OUT["subD"] * (152.2 + 7.122 * OUT["subT"]) # Effective heat capacity # Stability time step (CFL condition) - z_temp = OUT['subZ'][:, 1:] + z_temp = OUT["subZ"][:, 1:] c_eff_temp = c_eff[:, 1:] kk_temp = kk[:, 1:] dt_stab = ( - 0.5 * np.min(c_eff_temp, axis=1) * - np.min(z_temp, axis=1) ** 2 / - np.max(kk_temp, axis=1) / - C['dayseconds'] + 0.5 * np.min(c_eff_temp, axis=1) * np.min(z_temp, axis=1) ** 2 / np.max(kk_temp, axis=1) / C["dayseconds"] ) # ------ Heat Conduction Loop ------ - tt = np.zeros(grid['gpsum']) + tt = np.zeros(grid["gpsum"]) cond_dt_temp = np.zeros_like(tt, dtype=bool) - kdTdz = np.zeros_like(OUT['subT']) + kdTdz = np.zeros_like(OUT["subT"]) while np.any(tt < dt): - subT_old = OUT['subT'].copy() + subT_old = OUT["subT"].copy() dt_temp = np.minimum(dt_stab, dt - tt) tt += dt_temp cond_dt = dt_temp > 0 @@ -356,67 +359,53 @@ def heat_conduction(): # Calculate vertical heat fluxes kdTdz[cond_dt, 1] = ( - (kk[cond_dt, 0] * OUT['subZ'][cond_dt, 0] + - 0.5 * kk[cond_dt, 1] * OUT['subZ'][cond_dt, 1]) * - (subT_old[cond_dt, 1] - OUT['Tsurf'][cond_dt]) / - dz1[cond_dt] + (kk[cond_dt, 0] * OUT["subZ"][cond_dt, 0] + 0.5 * kk[cond_dt, 1] * OUT["subZ"][cond_dt, 1]) + * (subT_old[cond_dt, 1] - OUT["Tsurf"][cond_dt]) + / dz1[cond_dt] ) kdTdz[cond_dt, 2:] = ( - (kk[cond_dt, 1:-1] * OUT['subZ'][cond_dt, 1:-1] + - kk[cond_dt, 2:] * OUT['subZ'][cond_dt, 2:]) * - (subT_old[cond_dt, 2:] - subT_old[cond_dt, 1:-1]) / - dz2[cond_dt] + (kk[cond_dt, 1:-1] * OUT["subZ"][cond_dt, 1:-1] + kk[cond_dt, 2:] * OUT["subZ"][cond_dt, 2:]) + * (subT_old[cond_dt, 2:] - subT_old[cond_dt, 1:-1]) + / dz2[cond_dt] ) # Update layer-wise temperatures - C_day_dt = C['dayseconds'] * dt_temp[cond_dt] - OUT['subT'][cond_dt, 1] = ( - subT_old[cond_dt, 1] + - C_day_dt * - (kdTdz[cond_dt, 2] - kdTdz[cond_dt, 1]) / - (c_eff[cond_dt, 1] * - (0.5 * OUT['subZ'][cond_dt, 0] + - 0.5 * OUT['subZ'][cond_dt, 1] + - 0.25 * OUT['subZ'][cond_dt, 2])) + C_day_dt = C["dayseconds"] * dt_temp[cond_dt] + OUT["subT"][cond_dt, 1] = subT_old[cond_dt, 1] + C_day_dt * (kdTdz[cond_dt, 2] - kdTdz[cond_dt, 1]) / ( + c_eff[cond_dt, 1] + * (0.5 * OUT["subZ"][cond_dt, 0] + 0.5 * OUT["subZ"][cond_dt, 1] + 0.25 * OUT["subZ"][cond_dt, 2]) ) - OUT['subT'][cond_dt, 2:-1] = ( - subT_old[cond_dt, 2:-1] + - C_day_dt[:, np.newaxis] * - (kdTdz[cond_dt, 3:] - kdTdz[cond_dt, 2:-1]) / - ( - c_eff[cond_dt, 2:-1] * - (0.25 * OUT['subZ'][cond_dt, 1:-2] + - 0.5 * OUT['subZ'][cond_dt, 2:-1] + - 0.25 * OUT['subZ'][cond_dt, 3:]) - ) + OUT["subT"][cond_dt, 2:-1] = subT_old[cond_dt, 2:-1] + C_day_dt[:, np.newaxis] * ( + kdTdz[cond_dt, 3:] - kdTdz[cond_dt, 2:-1] + ) / ( + c_eff[cond_dt, 2:-1] + * ( + 0.25 * OUT["subZ"][cond_dt, 1:-2] + + 0.5 * OUT["subZ"][cond_dt, 2:-1] + + 0.25 * OUT["subZ"][cond_dt, 3:] + ) ) - OUT['subT'][cond_dt, -1] = ( - subT_old[cond_dt, -1] + - C_day_dt * - (C['geothermal_flux'] - kdTdz[cond_dt, -1]) / - ( - c_eff[cond_dt, -1] * - (0.25 * OUT['subZ'][cond_dt, -2] + - 0.75 * OUT['subZ'][cond_dt, -1]) - ) - ) + OUT["subT"][cond_dt, -1] = subT_old[cond_dt, -1] + C_day_dt * ( + C["geothermal_flux"] - kdTdz[cond_dt, -1] + ) / (c_eff[cond_dt, -1] * (0.25 * OUT["subZ"][cond_dt, -2] + 0.75 * OUT["subZ"][cond_dt, -1])) - OUT['subT'][:, 0] = ( - OUT['Tsurf'] + - (OUT['subT'][:, 1] - OUT['Tsurf']) / - (OUT['subZ'][:, 0] + 0.5 * OUT['subZ'][:, 1]) * - 0.5 * OUT['subZ'][:, 0] + OUT["subT"][:, 0] = ( + OUT["Tsurf"] + + (OUT["subT"][:, 1] - OUT["Tsurf"]) + / (OUT["subZ"][:, 0] + 0.5 * OUT["subZ"][:, 1]) + * 0.5 + * OUT["subZ"][:, 0] ) # Ensure temperatures do not exceed melting point - np.clip(OUT['subT'], None, C['T0'], out=OUT['subT']) + np.clip(OUT["subT"], None, C["T0"], out=OUT["subT"]) # Store effective conductivity and specific heat capacity - OUT['subCeff'] = c_eff - OUT['subK'] = kk + OUT["subCeff"] = c_eff + OUT["subK"] = kk return True @@ -424,66 +413,65 @@ def percolation_refreezing_and_storage(): ######################################################### # Percolation, refreezing and irreducible water storage ######################################################### - subW_old = OUT['subW'].copy() # Store the old water content + subW_old = OUT["subW"].copy() # Store the old water content gpsum, nl = OUT["subT"].shape # ------ Water Input ------ avail_W = ( - OUT['melt'] * 1e3 + # Meltwater - IN['rain'] * 1e3 + # Rainfall - (OUT['moist_condensation'] - OUT['moist_evaporation']) * 1e3 # Condensation or evaporation + OUT["melt"] * 1e3 # Meltwater + + IN["rain"] * 1e3 # Rainfall + + (OUT["moist_condensation"] - OUT["moist_evaporation"]) * 1e3 # Condensation or evaporation ) avail_W = np.maximum(avail_W, 0) # Ensure no negative water availability # ------ Refreezing and Irreducible Water Storage Limits ------ - OUT['cpi'] = 152.2 + 7.122 * OUT['subT'] # Specific heat capacity - c1 = OUT['cpi'] * OUT['subD'] * OUT['subZ'] * (C['T0'] - OUT['subT']) / C['Lm'] - c2 = OUT['subZ'] * (1 - OUT['subD'] / C['Dice']) * C['Dice'] + OUT["cpi"] = 152.2 + 7.122 * OUT["subT"] # Specific heat capacity + c1 = OUT["cpi"] * OUT["subD"] * OUT["subZ"] * (C["T0"] - OUT["subT"]) / C["Lm"] + c2 = OUT["subZ"] * (1 - OUT["subD"] / C["Dice"]) * C["Dice"] cond1 = c1 >= c2 # No need for separate cond2 as it's just the negation # Compute refreezing potential (`Wlim`) per layer Wlim = np.maximum(np.where(cond1, c2, c1), 0) # Maximum irreducible water storage (`mliqmax`) - mliqmax = np.zeros_like(OUT['subD']) - noice = OUT['subD'] < (C['Dice'] - 1) - factor = 3.3 * (C['Dice'] - OUT['subD'][noice]) / C['Dice'] + mliqmax = np.zeros_like(OUT["subD"]) + noice = OUT["subD"] < (C["Dice"] - 1) + factor = 3.3 * (C["Dice"] - OUT["subD"][noice]) / C["Dice"] exp_factor = np.exp(factor) irr_factor = 0.0143 * exp_factor / (1 - 0.0143 * exp_factor) mliqmax[noice] = ( - OUT['subD'][noice] * OUT['subZ'][noice] * - irr_factor * 0.05 * np.minimum(C['Dice'] - OUT['subD'][noice], 20) + OUT["subD"][noice] * OUT["subZ"][noice] * irr_factor * 0.05 * np.minimum(C["Dice"] - OUT["subD"][noice], 20) ) # Available irreducible water storage Wirr = mliqmax - subW_old # ------ Water Percolation ------ - z0 = C['perc_depth'] - zz = np.cumsum(OUT['subZ'], axis=1) - 0.5 * OUT['subZ'] - carrot = np.zeros_like(OUT['subZ']) + z0 = C["perc_depth"] + zz = np.cumsum(OUT["subZ"], axis=1) - 0.5 * OUT["subZ"] + carrot = np.zeros_like(OUT["subZ"]) - if phys['percolation'] == 'bucket': + if phys["percolation"] == "bucket": carrot[:, 0] = 1 # All water is added at the surface layer - elif phys['percolation'] == 'normal': - carrot = 2 * np.exp(-zz ** 2 / (2 * (z0 / 3) ** 2)) / (z0 / 3) / np.sqrt(2 * np.pi) - elif phys['percolation'] == 'linear': - carrot = 2 * (z0 - zz) / z0 ** 2 + elif phys["percolation"] == "normal": + carrot = 2 * np.exp(-(zz**2) / (2 * (z0 / 3) ** 2)) / (z0 / 3) / np.sqrt(2 * np.pi) + elif phys["percolation"] == "linear": + carrot = 2 * (z0 - zz) / z0**2 carrot = np.maximum(carrot, 0) - elif phys['percolation'] == 'uniform': + elif phys["percolation"] == "uniform": ind = np.argmin(np.abs(zz - z0), axis=1) - carrot[np.arange(carrot.shape[0]), :ind + 1] = 1 / z0 + carrot[np.arange(carrot.shape[0]), : ind + 1] = 1 / z0 else: raise ValueError("`phys['percolation']` is not set correctly!") - carrot *= OUT['subZ'] # Scale by layer thickness + carrot *= OUT["subZ"] # Scale by layer thickness carrot /= np.sum(carrot, axis=1)[:, np.newaxis] # Normalize per layer carrot *= avail_W[:, np.newaxis] # Distribute water input among layers # ------ Refreezing and Irreducible Water Storage Iteration ------ - RP = np.zeros_like(OUT['subZ']) # Refreezing potential - leftW = np.zeros(grid['gpsum']) # Remaining water - avail_W_loc = np.zeros(grid['gpsum']) # Available water per layer + RP = np.zeros_like(OUT["subZ"]) # Refreezing potential + leftW = np.zeros(grid["gpsum"]) # Remaining water + avail_W_loc = np.zeros(grid["gpsum"]) # Available water per layer for n in range(nl): # Compute available water per layer @@ -494,19 +482,19 @@ def percolation_refreezing_and_storage(): RP[cond1, n] = Wlim[cond1, n] # Refreezing limited by Wlim leftW[cond1] = avail_W_loc[cond1] - Wlim[cond1, n] - OUT['subW'][cond1, n] = subW_old[cond1, n] + np.minimum(leftW[cond1], Wirr[cond1, n]) + OUT["subW"][cond1, n] = subW_old[cond1, n] + np.minimum(leftW[cond1], Wirr[cond1, n]) # Condition: water <= refreezing limit RP[~cond1, n] = avail_W_loc[~cond1] - OUT['subW'][~cond1, n] = subW_old[~cond1, n] + OUT["subW"][~cond1, n] = subW_old[~cond1, n] # Update available water after refreezing - avail_W_loc -= RP[:, n] + (OUT['subW'][:, n] - subW_old[:, n]) + avail_W_loc -= RP[:, n] + (OUT["subW"][:, n] - subW_old[:, n]) # Update temperature and density after refreezing - refreeze_heat = C['Lm'] * RP[:, n] - OUT['subT'][:, n] += refreeze_heat / (OUT['subD'][:, n] * OUT['cpi'][:, n] * OUT['subZ'][:, n]) - OUT['subD'][:, n] += RP[:, n] / OUT['subZ'][:, n] + refreeze_heat = C["Lm"] * RP[:, n] + OUT["subT"][:, n] += refreeze_heat / (OUT["subD"][:, n] * OUT["cpi"][:, n] * OUT["subZ"][:, n]) + OUT["subD"][:, n] += RP[:, n] / OUT["subZ"][:, n] # Update leftover water avail_W = avail_W_loc @@ -517,8 +505,7 @@ def percolation_refreezing_and_storage(): # Calculate available pore space for storing slush water slushspace = np.maximum( - OUT["subZ"] * (1 - OUT["subD"] / C["Dice"]) * C["Dwater"] - OUT["subW"], - 0.0 + OUT["subZ"] * (1 - OUT["subD"] / C["Dice"]) * C["Dwater"] - OUT["subW"], 0.0 ) # shape: [grid.gpsum, nl] # Total available slush space across all layers @@ -526,18 +513,18 @@ def percolation_refreezing_and_storage(): # Calculate surface runoff (excess water at the surface) avail_W += np.sum(OUT["subS"], axis=1) # Update available water with slush from all layers - OUT['runoff_surface'] = np.maximum(avail_W - total_slushspace, 0.0) # Excess water at the surface + OUT["runoff_surface"] = np.maximum(avail_W - total_slushspace, 0.0) # Excess water at the surface # Update slush water content after new water input and runoff avail_S = np.minimum(avail_W, total_slushspace) # Available slush water for storage # Calculate slush water runoff and reduce avail_S accordingly - OUT['runoff_slush'] = avail_S - 1.0 / (1.0 + dt / C["Trunoff"]) * avail_S + OUT["runoff_slush"] = avail_S - 1.0 / (1.0 + dt / C["Trunoff"]) * avail_S avail_S = 1.0 / (1.0 + dt / C["Trunoff"]) * avail_S avail_S[avail_S < 1e-25] = 0.0 # Set near-zero available slush to zero # Initialize slush water in all layers - OUT["subS"] = np.zeros((grid['gpsum'], nl)) + OUT["subS"] = np.zeros((grid["gpsum"], nl)) # Bottom-up filling of pore space with slush water for n in range(nl - 1, -1, -1): # Loop from bottom (nl) to top (1) @@ -551,15 +538,15 @@ def percolation_refreezing_and_storage(): ##################################### # Determine whether cold content or density limits the amount of refreezing OUT["cpi"] = 152.2 + 7.122 * OUT["subT"] - c1 = OUT["cpi"] * OUT["subD"] * OUT["subZ"] * (C['T0'] - OUT["subT"]) / C['Lm'] - c2 = OUT["subZ"] * (1 - OUT["subD"] / C['Dice']) * C['Dice'] + c1 = OUT["cpi"] * OUT["subD"] * OUT["subZ"] * (C["T0"] - OUT["subT"]) / C["Lm"] + c2 = OUT["subZ"] * (1 - OUT["subD"] / C["Dice"]) * C["Dice"] Wlim = np.minimum(c1, c2) # Available slush water slush_W = OUT["subS"].copy() # Make a copy to avoid overwriting inputs (vectorized across all layers) # Determine refreezing amounts - layer_cond = (OUT["subS"] > 0) & (OUT["subT"] < C['T0']) + layer_cond = (OUT["subS"] > 0) & (OUT["subT"] < C["T0"]) Wlim_effective = Wlim * layer_cond slush_W_effective = slush_W * layer_cond RS = np.where(slush_W_effective > Wlim_effective, Wlim_effective, slush_W_effective) @@ -568,7 +555,7 @@ def percolation_refreezing_and_storage(): OUT["subS"] -= RS # Update temperature after refreezing - OUT["subT"] += (C['Lm'] * RS) / (OUT["subD"] * OUT["cpi"] * OUT["subZ"]) + OUT["subT"] += (C["Lm"] * RS) / (OUT["subD"] * OUT["cpi"] * OUT["subZ"]) # Update density after refreezing OUT["subD"] += RS / OUT["subZ"] @@ -577,27 +564,27 @@ def percolation_refreezing_and_storage(): ### REFREEZING OF IRREDUCIBLE WATER ############################################################# # Determine whether cold content or density limits the amount of refreezing - OUT['cpi'] = 152.2 + 7.122 * OUT['subT'] - c1 = OUT['cpi'] * OUT['subD'] * OUT['subZ'] * (C['T0'] - OUT['subT']) / C['Lm'] - c2 = OUT['subZ'] * (1 - OUT['subD'] / C['Dice']) * C['Dice'] + OUT["cpi"] = 152.2 + 7.122 * OUT["subT"] + c1 = OUT["cpi"] * OUT["subD"] * OUT["subZ"] * (C["T0"] - OUT["subT"]) / C["Lm"] + c2 = OUT["subZ"] * (1 - OUT["subD"] / C["Dice"]) * C["Dice"] Wlim = np.where(c1 >= c2, c2, c1) # Calculate refreezing amounts - valid_mask = (OUT['subW'] > 0) & (OUT['subT'] < C['T0']) - RI = np.minimum(OUT['subW'], Wlim) * valid_mask + valid_mask = (OUT["subW"] > 0) & (OUT["subT"] < C["T0"]) + RI = np.minimum(OUT["subW"], Wlim) * valid_mask # Update water content (subW), temperature (subT), and density (subD) - OUT['subW'] -= RI - OUT['subT'] += (C['Lm'] * RI) / (OUT['subD'] * OUT['cpi'] * OUT['subZ']) - OUT['subD'] += RI / OUT['subZ'] + OUT["subW"] -= RI + OUT["subT"] += (C["Lm"] * RI) / (OUT["subD"] * OUT["cpi"] * OUT["subZ"]) + OUT["subD"] += RI / OUT["subZ"] # Determine total refreezing and individual contributions, and total slush water and irreducible water - OUT['refr'] = 1e-3 * (np.sum(RP, axis=1) + np.sum(RS, axis=1) + np.sum(RI, axis=1)) - OUT['refr_P'] = 1e-3 * np.sum(RP, axis=1) # Refreezing of percolating water - OUT['refr_S'] = 1e-3 * np.sum(RS, axis=1) # Refreezing of slush water - OUT['refr_I'] = 1e-3 * np.sum(RI, axis=1) # Refreezing of irreducible water - OUT['slushw'] = np.sum(OUT['subS'], axis=1) # Total stored slush water - OUT['irrw'] = np.sum(OUT['subW'], axis=1) # Total stored irreducible water + OUT["refr"] = 1e-3 * (np.sum(RP, axis=1) + np.sum(RS, axis=1) + np.sum(RI, axis=1)) + OUT["refr_P"] = 1e-3 * np.sum(RP, axis=1) # Refreezing of percolating water + OUT["refr_S"] = 1e-3 * np.sum(RS, axis=1) # Refreezing of slush water + OUT["refr_I"] = 1e-3 * np.sum(RI, axis=1) # Refreezing of irreducible water + OUT["slushw"] = np.sum(OUT["subS"], axis=1) # Total stored slush water + OUT["irrw"] = np.sum(OUT["subW"], axis=1) # Total stored irreducible water return True @@ -611,45 +598,41 @@ def layer_merging_and_splitting(): subW_old = OUT["subW"].copy() subT_old = OUT["subT"].copy() subS_old = OUT["subS"].copy() - for n in range(len(grid['split'])): # Iterate through split points - split = grid['split'][n] - 1 + for n in range(len(grid["split"])): # Iterate through split points + split = grid["split"][n] - 1 # Merge Layers (Accumulation Case) - cond_merge = (OUT["subZ"][:, split] <= (2.0 ** n) * grid['max_subZ']) & \ - (grid['mask'] == 1) + cond_merge = (OUT["subZ"][:, split] <= (2.0**n) * grid["max_subZ"]) & (grid["mask"] == 1) # Update merged layers OUT["subZ"][cond_merge, split - 1] = subZ_old[cond_merge, split - 1] + subZ_old[cond_merge, split] OUT["subW"][cond_merge, split - 1] = subW_old[cond_merge, split - 1] + subW_old[cond_merge, split] OUT["subS"][cond_merge, split - 1] = subS_old[cond_merge, split - 1] + subS_old[cond_merge, split] OUT["subD"][cond_merge, split - 1] = ( - (subZ_old[cond_merge, split - 1] * subD_old[cond_merge, split - 1] + - subZ_old[cond_merge, split] * subD_old[cond_merge, split]) / - (subZ_old[cond_merge, split - 1] + subZ_old[cond_merge, split]) - ) + subZ_old[cond_merge, split - 1] * subD_old[cond_merge, split - 1] + + subZ_old[cond_merge, split] * subD_old[cond_merge, split] + ) / (subZ_old[cond_merge, split - 1] + subZ_old[cond_merge, split]) OUT["subT"][cond_merge, split - 1] = ( - (subZ_old[cond_merge, split - 1] * subT_old[cond_merge, split - 1] + - subZ_old[cond_merge, split] * subT_old[cond_merge, split]) / - (subZ_old[cond_merge, split - 1] + subZ_old[cond_merge, split]) - ) + subZ_old[cond_merge, split - 1] * subT_old[cond_merge, split - 1] + + subZ_old[cond_merge, split] * subT_old[cond_merge, split] + ) / (subZ_old[cond_merge, split - 1] + subZ_old[cond_merge, split]) # Shift properties up for merged layers - OUT["subZ"][cond_merge, split:-1] = subZ_old[cond_merge, split + 1:] - OUT["subW"][cond_merge, split:-1] = subW_old[cond_merge, split + 1:] - OUT["subS"][cond_merge, split:-1] = subS_old[cond_merge, split + 1:] - OUT["subD"][cond_merge, split:-1] = subD_old[cond_merge, split + 1:] - OUT["subT"][cond_merge, split:-1] = subT_old[cond_merge, split + 1:] + OUT["subZ"][cond_merge, split:-1] = subZ_old[cond_merge, split + 1 :] + OUT["subW"][cond_merge, split:-1] = subW_old[cond_merge, split + 1 :] + OUT["subS"][cond_merge, split:-1] = subS_old[cond_merge, split + 1 :] + OUT["subD"][cond_merge, split:-1] = subD_old[cond_merge, split + 1 :] + OUT["subT"][cond_merge, split:-1] = subT_old[cond_merge, split + 1 :] # Adjust the newly added layer at the top - OUT["subZ"][cond_merge, -1] = 2.0 ** len(grid['split']) * grid['max_subZ'] + OUT["subZ"][cond_merge, -1] = 2.0 ** len(grid["split"]) * grid["max_subZ"] OUT["subT"][cond_merge, -1] = 2.0 * subT_old[cond_merge, -1] - subT_old[cond_merge, -2] OUT["subD"][cond_merge, -1] = subD_old[cond_merge, -1] OUT["subW"][cond_merge, -1] = 0.0 OUT["subS"][cond_merge, -1] = 0.0 # Split Layers (Ablation Case) - cond_split = (OUT["subZ"][:, split - 2] > (2.0 ** n) * grid['max_subZ']) & \ - (grid['mask'] == 1) + cond_split = (OUT["subZ"][:, split - 2] > (2.0**n) * grid["max_subZ"]) & (grid["mask"] == 1) # Update split layers OUT["subZ"][cond_split, split - 2] *= 0.5 @@ -665,15 +648,15 @@ def layer_merging_and_splitting(): OUT["subD"][cond_split, split - 1] = OUT["subD"][cond_split, split - 2] # Shift properties down for split layers - OUT["subZ"][cond_split, split:-1] = subZ_old[cond_split, split - 1:-2] - OUT["subW"][cond_split, split:-1] = subW_old[cond_split, split - 1:-2] - OUT["subS"][cond_split, split:-1] = subS_old[cond_split, split - 1:-2] - OUT["subT"][cond_split, split:-1] = subT_old[cond_split, split - 1:-2] - OUT["subD"][cond_split, split:-1] = subD_old[cond_split, split - 1:-2] + OUT["subZ"][cond_split, split:-1] = subZ_old[cond_split, split - 1 : -2] + OUT["subW"][cond_split, split:-1] = subW_old[cond_split, split - 1 : -2] + OUT["subS"][cond_split, split:-1] = subS_old[cond_split, split - 1 : -2] + OUT["subT"][cond_split, split:-1] = subT_old[cond_split, split - 1 : -2] + OUT["subD"][cond_split, split:-1] = subD_old[cond_split, split - 1 : -2] # Update runoff contributions - OUT['runoff_irr_deep'][cond_split] += subW_old[cond_split, -1] - OUT['runoff_slush'][cond_split] += subS_old[cond_split, -1] + OUT["runoff_irr_deep"][cond_split] += subW_old[cond_split, -1] + OUT["runoff_slush"][cond_split] += subS_old[cond_split, -1] return True @@ -683,12 +666,14 @@ def runoff(): ########################################### # Update runoff of irreducible water below the model bottom - OUT["runoff_irr_deep_mean"] = OUT["runoff_irr_deep_mean"] * (1 - dt / C["yeardays"]) + \ - OUT['runoff_irr_deep'] * (dt / C["yeardays"]) + OUT["runoff_irr_deep_mean"] = OUT["runoff_irr_deep_mean"] * (1 - dt / C["yeardays"]) + OUT[ + "runoff_irr_deep" + ] * (dt / C["yeardays"]) # Calculate total runoff [in meters water equivalent per timestep] OUT["runoff"] = 1e-3 * ( - OUT["runoff_surface"] + OUT["runoff_slush"] + OUT["runoff_irr"] + OUT["runoff_irr_deep_mean"]) + OUT["runoff_surface"] + OUT["runoff_slush"] + OUT["runoff_irr"] + OUT["runoff_irr_deep_mean"] + ) # Surface runoff [in meters water equivalent per timestep] OUT["runoff_surf"] = 1e-3 * OUT["runoff_surface"] @@ -697,7 +682,7 @@ def runoff(): OUT["runoff_slush"] = 1e-3 * OUT["runoff_slush"] # Irreducible water runoff within the domain [in meters water equivalent per timestep] - OUT["runoff_irr"] = 1e-3 * OUT['runoff_irr'] + OUT["runoff_irr"] = 1e-3 * OUT["runoff_irr"] # Irreducible water runoff below the base of the domain [in meters water equivalent per timestep] OUT["runoff_irr_deep"] = 1e-3 * OUT["runoff_irr_deep_mean"] @@ -711,6 +696,6 @@ def runoff(): percolation_refreezing_and_storage() layer_merging_and_splitting() runoff() - OUT['T_ice'] = OUT['subT'][:, -1] + OUT["T_ice"] = OUT["subT"][:, -1] return OUT diff --git a/src/ebfm/LOOP_climate_forcing.py b/src/ebfm/LOOP_climate_forcing.py index 8552c3d2..82df1283 100644 --- a/src/ebfm/LOOP_climate_forcing.py +++ b/src/ebfm/LOOP_climate_forcing.py @@ -2,6 +2,7 @@ from coupling import Coupler + def main(C, grid, IN, t, time, OUT, cpl: Coupler) -> tuple[dict, dict]: """ Meteorological forcing: Specify or read meteorological input and derive @@ -38,53 +39,52 @@ def main(C, grid, IN, t, time, OUT, cpl: Coupler) -> tuple[dict, dict]: ########################################################### # Annual snow accumulation - OUT['ys'] = (1.0 - (1.0 / (C['yeardays'] / time['dt']))) * OUT['ys'] + IN['P'] * 1e3 - logys = np.log(OUT['ys']) - IN['yearsnow'] = np.tile(OUT['ys'][:, np.newaxis], (1, grid['nl'])) - IN['logyearsnow'] = np.tile(logys[:, np.newaxis], (1, grid['nl'])) + OUT["ys"] = (1.0 - (1.0 / (C["yeardays"] / time["dt"]))) * OUT["ys"] + IN["P"] * 1e3 + logys = np.log(OUT["ys"]) + IN["yearsnow"] = np.tile(OUT["ys"][:, np.newaxis], (1, grid["nl"])) + IN["logyearsnow"] = np.tile(logys[:, np.newaxis], (1, grid["nl"])) # Vapor pressure, relative and specific humidity - VPsat = ( - C['VP0'] * np.exp(C['Lv'] / C['Rv'] * (1.0 / 273.15 - 1.0 / IN['T'])) * (IN['T'] >= 273.15) + - C['VP0'] * np.exp(C['Ls'] / C['Rv'] * (1.0 / 273.15 - 1.0 / IN['T'])) * (IN['T'] < 273.15) - ) - - if cpl.couple_to_icon_atmo: # q from ICON, calculate VP and RH - IN['VP'] = IN['q'] * IN['Pres'] / C['eps'] - IN['RH'] = IN['VP'] / VPsat - else: # RH from input, calculate VP and q - IN['VP'] = IN['RH'] * VPsat - IN['q'] = IN['RH'] * (VPsat * C['eps'] / IN['Pres']) + VPsat = C["VP0"] * np.exp(C["Lv"] / C["Rv"] * (1.0 / 273.15 - 1.0 / IN["T"])) * (IN["T"] >= 273.15) + C[ + "VP0" + ] * np.exp(C["Ls"] / C["Rv"] * (1.0 / 273.15 - 1.0 / IN["T"])) * (IN["T"] < 273.15) + + if cpl.couple_to_icon_atmo: # q from ICON, calculate VP and RH + IN["VP"] = IN["q"] * IN["Pres"] / C["eps"] + IN["RH"] = IN["VP"] / VPsat + else: # RH from input, calculate VP and q + IN["VP"] = IN["RH"] * VPsat + IN["q"] = IN["RH"] * (VPsat * C["eps"] / IN["Pres"]) # Air density - IN['Dair'] = IN['Pres'] / (C['Rd'] * IN['T']) + IN["Dair"] = IN["Pres"] / (C["Rd"] * IN["T"]) # Time since last snowfall event - snowfall_mask = (IN['snow'] / (time['dt'] * 24 * 3600)) > C['Pthres'] - OUT['timelastsnow'][snowfall_mask] = time['TCUR'] + snowfall_mask = (IN["snow"] / (time["dt"] * 24 * 3600)) > C["Pthres"] + OUT["timelastsnow"][snowfall_mask] = time["TCUR"] if t == 1: - OUT['timelastsnow'][:] = time['TCUR'] + OUT["timelastsnow"][:] = time["TCUR"] # Potential temperature and lapse rate - IN['Theta'] = IN['T'] * (C['Pref'] / IN['Pres']) ** (C['Rd'] / C['Cp']) - all_same = np.all(grid['z'] == grid['z'][0]) + IN["Theta"] = IN["T"] * (C["Pref"] / IN["Pres"]) ** (C["Rd"] / C["Cp"]) + all_same = np.all(grid["z"] == grid["z"][0]) if not all_same: - poly_coeff = np.polyfit(grid['z'], IN['Theta'], deg=1) - IN['Theta_lapse'] = max(poly_coeff[0], 0.0015) + poly_coeff = np.polyfit(grid["z"], IN["Theta"], deg=1) + IN["Theta_lapse"] = max(poly_coeff[0], 0.0015) else: - IN['Theta_lapse'] = 0.0015 + IN["Theta_lapse"] = 0.0015 ########################################################### # STORE RELEVANT VARIABLES IN OUT ########################################################### - OUT['climT'] = IN['T'] - OUT['climP'] = IN['P'] - OUT['climC'] = IN['C'] - OUT['climRH'] = IN['RH'] - OUT['climWS'] = IN['WS'] - OUT['climPres'] = IN['Pres'] - OUT['climsnow'] = IN['snow'] - OUT['climrain'] = IN['rain'] + OUT["climT"] = IN["T"] + OUT["climP"] = IN["P"] + OUT["climC"] = IN["C"] + OUT["climRH"] = IN["RH"] + OUT["climWS"] = IN["WS"] + OUT["climPres"] = IN["Pres"] + OUT["climsnow"] = IN["snow"] + OUT["climrain"] = IN["rain"] return IN, OUT @@ -105,44 +105,42 @@ def set_random_weather_data(IN, C, time, grid): ############################## # Example: Random Conditions ############################## - yearfrac = time['TCUR'].timetuple().tm_yday / C['yeardays'] + yearfrac = time["TCUR"].timetuple().tm_yday / C["yeardays"] # Air temperature (K) T_amplitude = 10.0 # Seasonal temperature amplitude (K) T_mean_sea_level = 268.0 # Mean sea level temperature (K) T_lapse_rate = -0.005 # Temperature lapse rate (K m-1) - IN['T'] = ( - T_mean_sea_level + T_amplitude * np.sin(2 * np.pi * yearfrac - 0.65 * np.pi) - ) - IN['T'] += T_lapse_rate * grid['z'] + IN["T"] = T_mean_sea_level + T_amplitude * np.sin(2 * np.pi * yearfrac - 0.65 * np.pi) + IN["T"] += T_lapse_rate * grid["z"] # Precipitation (m w.e.) P_annual_sea_level = 0.5 # Annual precipitation at sea level (m w.e.) P_z_gradient = 0.1 # Precipitation - elevation gradient (% m-1) - if time['TCUR'].isocalendar()[1] != time['TPREV'].isocalendar()[1]: - IN['P'][:] = (P_annual_sea_level / 52.0) * (1 + P_z_gradient * grid['z'] / 100.0) + if time["TCUR"].isocalendar()[1] != time["TPREV"].isocalendar()[1]: + IN["P"][:] = (P_annual_sea_level / 52.0) * (1 + P_z_gradient * grid["z"] / 100.0) else: - IN['P'][:] = 0.0 + IN["P"][:] = 0.0 # Cloud cover (fraction) - IN['C'][:] = 1.0 if time['TCUR'].isocalendar()[1] % 2 == 0 else 0.0 + IN["C"][:] = 1.0 if time["TCUR"].isocalendar()[1] % 2 == 0 else 0.0 # Relative humidity (fraction) - IN['RH'][:] = 0.8 if time['TCUR'].isocalendar()[1] % 2 == 0 else 0.5 + IN["RH"][:] = 0.8 if time["TCUR"].isocalendar()[1] % 2 == 0 else 0.5 # Wind speed (m s-1) max_WS = 10.0 # Max wind speed - IN['WS'][:] = np.random.uniform(0.0, max_WS, size=grid['gpsum']) + IN["WS"][:] = np.random.uniform(0.0, max_WS, size=grid["gpsum"]) # Air pressure (Pa) Pres_sea_level = 1015e2 # Sea level pressure (Pa) - IN['Pres'][:] = Pres_sea_level * np.exp(-1.244e-4 * grid['z']) + IN["Pres"][:] = Pres_sea_level * np.exp(-1.244e-4 * grid["z"]) # Snowfall and rainfall - IN['snow'] = IN['P'] * (IN['T'] < C['rainsnowT'] - 1) - IN['rain'] = IN['P'] * (IN['T'] > C['rainsnowT'] + 1) - in_between_mask = (IN['T'] < C['rainsnowT'] + 1) & (IN['T'] > C['rainsnowT'] - 1) - IN['snow'] += IN['P'] * (C['rainsnowT'] - IN['T'] + 1) / 2 * in_between_mask - IN['rain'] += IN['P'] * (1 + IN['T'] - C['rainsnowT']) / 2 * in_between_mask + IN["snow"] = IN["P"] * (IN["T"] < C["rainsnowT"] - 1) + IN["rain"] = IN["P"] * (IN["T"] > C["rainsnowT"] + 1) + in_between_mask = (IN["T"] < C["rainsnowT"] + 1) & (IN["T"] > C["rainsnowT"] - 1) + IN["snow"] += IN["P"] * (C["rainsnowT"] - IN["T"] + 1) / 2 * in_between_mask + IN["rain"] += IN["P"] * (1 + IN["T"] - C["rainsnowT"]) / 2 * in_between_mask return IN diff --git a/src/ebfm/LOOP_general_functions.py b/src/ebfm/LOOP_general_functions.py index 2c0f6c71..e4dec768 100644 --- a/src/ebfm/LOOP_general_functions.py +++ b/src/ebfm/LOOP_general_functions.py @@ -19,14 +19,15 @@ def print_time(t, time): # time['TE'] = datetime.strptime(time['te'], '%Y-%m-%d %H:%M:%S') # End time # Calculate current and previous model times based on the time step - time['TCUR'] = time['ts'] + timedelta(days=(t - 1) * time['dt']) - time['TPREV'] = time['ts'] + timedelta(days=(t - 2) * time['dt']) + time["TCUR"] = time["ts"] + timedelta(days=(t - 1) * time["dt"]) + time["TPREV"] = time["ts"] + timedelta(days=(t - 2) * time["dt"]) # Extract year, month, day, and hour from the current time - tempdate = time['TCUR'].timetuple() + tempdate = time["TCUR"].timetuple() # Print the current time parameters print( - f'Year: {tempdate.tm_year:4d} Month: {tempdate.tm_mon:2d} Day: {tempdate.tm_mday:2d} Hour: {tempdate.tm_hour:2d}') + f"Year: {tempdate.tm_year:4d} Month: {tempdate.tm_mon:2d} Day: {tempdate.tm_mday:2d} Hour: {tempdate.tm_hour:2d}" + ) return time diff --git a/src/ebfm/LOOP_mass_balance.py b/src/ebfm/LOOP_mass_balance.py index 9781306b..2afb5233 100644 --- a/src/ebfm/LOOP_mass_balance.py +++ b/src/ebfm/LOOP_mass_balance.py @@ -15,9 +15,15 @@ def main(OUT, IN, C): """ # Climatic mass balance - OUT["smb"] = (IN["snow"] + IN["rain"] - OUT["runoff"] - + OUT["moist_deposition"] + OUT["moist_condensation"] - - OUT["moist_sublimation"] - OUT["moist_evaporation"]) + OUT["smb"] = ( + IN["snow"] + + IN["rain"] + - OUT["runoff"] + + OUT["moist_deposition"] + + OUT["moist_condensation"] + - OUT["moist_sublimation"] + - OUT["moist_evaporation"] + ) OUT["smb_cumulative"] += OUT["smb"] diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index d9fcce0a..a30ca699 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -39,11 +39,12 @@ def main(OUTFILE, io, OUT, grid, t, time, C): ["subT", "K", "sample", "Temperature"], ["subS", "mm w.e.", "sample", "Slush water content"], ["subW", "mm w.e.", "sample", "Irreducible water"], - ["subZ", "m", "sample", "Layer thickness"] + ["subZ", "m", "sample", "Layer thickness"], ] - io["varsout"] = [{"varname": v[0], "units": v[1], "type": v[2], "description": v[3]} for v in - OUTFILE["varsout"]] + io["varsout"] = [ + {"varname": v[0], "units": v[1], "type": v[2], "description": v[3]} for v in OUTFILE["varsout"] + ] # Update OUTFILE.TEMP with variables to be stored for entry in OUTFILE["varsout"]: @@ -64,7 +65,6 @@ def main(OUTFILE, io, OUT, grid, t, time, C): elif var_type == "sum": OUTFILE["TEMP"][varname] += temp_long - def save_binary_files(): """ Write model output to binary files and save run information. @@ -106,12 +106,7 @@ def save_binary_files(): file.close() # Prepare the runinfo dictionary - runinfo = { - "grid": grid, - "time": time, - "IOout": io, - "Cout": C - } + runinfo = {"grid": grid, "time": time, "IOout": io, "Cout": C} ################ WORK IN PROGRESS ############### @@ -186,9 +181,7 @@ def save_netCDF_file(): nc_var.description = var_desc # Define a time variable to track simulation steps - nc_time = io["nc_file"].createVariable( - "time", np.float64, ("time",), zlib=True, fill_value=-9999.0 - ) + nc_time = io["nc_file"].createVariable("time", np.float64, ("time",), zlib=True, fill_value=-9999.0) nc_time.units = time_units nc_time.calendar = time_calendar nc_time.description = "Time at which data is recorded, in days since 1970-01-01 00:00:00" @@ -212,9 +205,7 @@ def save_netCDF_file(): if varname.startswith("sub"): var_3D = np.full((grid["x_2D"].size, grid["nl"]), -9999.0) var_3D[grid["ind"], :] = var_1D - var_4D = var_3D.reshape(-1, grid["nl"]).reshape( - *grid["x_2D"].shape, grid["nl"] - ) + var_4D = var_3D.reshape(-1, grid["nl"]).reshape(*grid["x_2D"].shape, grid["nl"]) io["nc_file"][varname][time_index, :, :, :] = var_4D else: var_2D = np.full(grid["x_2D"].shape, -9999.0) diff --git a/src/ebfm/grid.py b/src/ebfm/grid.py index ae770426..8201f594 100644 --- a/src/ebfm/grid.py +++ b/src/ebfm/grid.py @@ -1,6 +1,7 @@ from enum import Enum + class GridInputType(Enum): - MATLAB = 'matlab' # .mat grid file with elevation - CUSTOM = 'custom' # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file - ELMER = 'elmer' # Elmer/Ice mesh file with elevation in z-coordinate \ No newline at end of file + MATLAB = "matlab" # .mat grid file with elevation + CUSTOM = "custom" # Elmer/Ice mesh file for xy-coordinates and separate NetCDF elevation file + ELMER = "elmer" # Elmer/Ice mesh file with elevation in z-coordinate diff --git a/src/elmer/mesh.py b/src/elmer/mesh.py index c4cfe12e..b67e91a4 100644 --- a/src/elmer/mesh.py +++ b/src/elmer/mesh.py @@ -2,21 +2,35 @@ from numpy.typing import NDArray import pyproj -class Mesh(): + +class Mesh: """A generic 3D Mesh""" - x_vertices: NDArray[np.float64] # x-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] - y_vertices: NDArray[np.float64] # y-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] - z_vertices: NDArray[np.float64] # z-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] - lon: NDArray[np.float64] # longitude coordinates of vertices in radians, ordering follows local ids [0,1,...,n_vertices-1] - lat: NDArray[np.float64] # latitude coordinates of vertices in radians, ordering follows local ids [0,1,...,n_vertices-1] + + x_vertices: NDArray[np.float64] # x-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] + y_vertices: NDArray[np.float64] # y-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] + z_vertices: NDArray[np.float64] # z-coordinates of vertices, ordering follows local ids [0,1,...,n_vertices-1] + lon: NDArray[ + np.float64 + ] # longitude coordinates of vertices in radians, ordering follows local ids [0,1,...,n_vertices-1] + lat: NDArray[ + np.float64 + ] # latitude coordinates of vertices in radians, ordering follows local ids [0,1,...,n_vertices-1] # @TODO later add slope # dzdx: NDArray[np.float64] # z-slope in x-direction # dzdy: NDArray[np.float64] # z-slope in y-direction - vertex_ids: NDArray[np.int_] # IDs of vertices, ordering follows local ids [0,1,...,n_vertices-1] + vertex_ids: NDArray[np.int_] # IDs of vertices, ordering follows local ids [0,1,...,n_vertices-1] cell_to_vertex: NDArray[np.int_] # Mapping from cells to their local vertex IDs - cell_ids: NDArray[np.int_] # IDs of cells (triangles), ordering follows local ids [0,1,...,n_cells-1] + cell_ids: NDArray[np.int_] # IDs of cells (triangles), ordering follows local ids [0,1,...,n_cells-1] - def __init__(self, x_vertices: NDArray[np.float64], y_vertices: NDArray[np.float64], z_vertices: NDArray[np.float64], cell_to_vertex: NDArray[np.int_], vertex_ids: NDArray[np.int_], cell_ids: NDArray[np.int_]): + def __init__( + self, + x_vertices: NDArray[np.float64], + y_vertices: NDArray[np.float64], + z_vertices: NDArray[np.float64], + cell_to_vertex: NDArray[np.int_], + vertex_ids: NDArray[np.int_], + cell_ids: NDArray[np.int_], + ): self.x_vertices = x_vertices self.y_vertices = y_vertices self.z_vertices = z_vertices @@ -27,10 +41,20 @@ def __init__(self, x_vertices: NDArray[np.float64], y_vertices: NDArray[np.float transformer = pyproj.Transformer.from_crs(3413, 4326, always_xy=True) self.lon, self.lat = transformer.transform(self.x_vertices, self.y_vertices, radians=True) + class TriangleMesh(Mesh): """A 3D Mesh consisting of triangular elements.""" + num_vertices_per_cell = 3 - def __init__(self, x_vertices: NDArray[np.float64], y_vertices: NDArray[np.float64], z_vertices: NDArray[np.float64], cell_to_vertex: NDArray[np.int_], vertex_ids: NDArray[np.int_], cell_ids: NDArray[np.int_]): - assert(cell_to_vertex.shape[1] == self.num_vertices_per_cell) # a triangle mesh has 3 nodes for all cells - super(TriangleMesh, self).__init__(x_vertices, y_vertices, z_vertices, cell_to_vertex, vertex_ids, cell_ids) \ No newline at end of file + def __init__( + self, + x_vertices: NDArray[np.float64], + y_vertices: NDArray[np.float64], + z_vertices: NDArray[np.float64], + cell_to_vertex: NDArray[np.int_], + vertex_ids: NDArray[np.int_], + cell_ids: NDArray[np.int_], + ): + assert cell_to_vertex.shape[1] == self.num_vertices_per_cell # a triangle mesh has 3 nodes for all cells + super(TriangleMesh, self).__init__(x_vertices, y_vertices, z_vertices, cell_to_vertex, vertex_ids, cell_ids) diff --git a/src/elmer/parser.py b/src/elmer/parser.py index 834fd93e..163e835d 100644 --- a/src/elmer/parser.py +++ b/src/elmer/parser.py @@ -7,6 +7,7 @@ """This file contains parsers for different file formats used in Elmer """ + def parse_header(header_file: Path) -> tuple[int, int]: """Parse Elmer mesh header file. Args: @@ -22,23 +23,32 @@ def parse_header(header_file: Path) -> tuple[int, int]: # 3. Number of 202 boundary edges (will be ignored) n_expected_columns: int = 3 if df.shape[1] != n_expected_columns: - raise ValueError(f"The header file has only {df.shape[1]} columns, but it should have exactly {n_expected_columns}.") - + raise ValueError( + f"The header file has only {df.shape[1]} columns, but it should have exactly {n_expected_columns}." + ) + if isinstance(df.iloc[0, 0], np.int_): n_vertices: int = int(df.iloc[0, 0]) else: - raise ValueError(f"Expected the first column of the header file to contain an integer, but got {type(df.iloc[0, 0])}.") - + raise ValueError( + f"Expected the first column of the header file to contain an integer, but got {type(df.iloc[0, 0])}." + ) + if isinstance(df.iloc[0, 1], np.int_): n_cells: int = int(df.iloc[0, 1]) else: - raise ValueError(f"Expected the second column of the header file to contain an integer, but got {type(df.iloc[0, 1])}.") + raise ValueError( + f"Expected the second column of the header file to contain an integer, but got {type(df.iloc[0, 1])}." + ) # Only the first row of the header file is relevant for the counts; other rows will be ignored. - + return n_vertices, n_cells -def parse_nodes(nodes_file: Path) -> tuple[NDArray[np.int_], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: + +def parse_nodes( + nodes_file: Path, +) -> tuple[NDArray[np.int_], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: """Parse Elmer mesh nodes file. Args: nodes_file (Path): Path to the Elmer mesh nodes file. @@ -51,19 +61,22 @@ def parse_nodes(nodes_file: Path) -> tuple[NDArray[np.int_], NDArray[np.float64] # 1. Global Node ID # 2. Node type (usually 0 for regular nodes) # 3. x coordinate - # 4. y coordinate + # 4. y coordinate # 5. z coordinate n_expected_columns: int = 5 if df.shape[1] < n_expected_columns: - raise ValueError(f"The nodes file has only {df.shape[1]} columns, but it should have exactly {n_expected_columns}.") + raise ValueError( + f"The nodes file has only {df.shape[1]} columns, but it should have exactly {n_expected_columns}." + ) - global_ids : NDArray[np.int_] = df.iloc[:, 0].to_numpy() - x : NDArray[np.float64] = df.iloc[:, 2].to_numpy() - y : NDArray[np.float64] = df.iloc[:, 3].to_numpy() - z : NDArray[np.float64] = df.iloc[:, 4].to_numpy() + global_ids: NDArray[np.int_] = df.iloc[:, 0].to_numpy() + x: NDArray[np.float64] = df.iloc[:, 2].to_numpy() + y: NDArray[np.float64] = df.iloc[:, 3].to_numpy() + z: NDArray[np.float64] = df.iloc[:, 4].to_numpy() return global_ids, x, y, z + def parse_elements(elements_file: Path) -> tuple[NDArray[np.int_], NDArray[np.int_]]: """Parse Elmer mesh elements file. Args: @@ -71,7 +84,7 @@ def parse_elements(elements_file: Path) -> tuple[NDArray[np.int_], NDArray[np.in Returns: tuple: Element IDs and connectivity. """ - df = pd.read_csv(elements_file, sep="\\s+", header=None, dtype={0:int, 3:float, 4:float, 5:float}) + df = pd.read_csv(elements_file, sep="\\s+", header=None, dtype={0: int, 3: float, 4: float, 5: float}) # A valid mesh.elements file in Elmer should have at least 2 columns: # 1. Global Element ID @@ -79,9 +92,11 @@ def parse_elements(elements_file: Path) -> tuple[NDArray[np.int_], NDArray[np.in # 4. - 6. column connectivity (list of node IDs that form the triangular element) n_expected_columns: int = 6 if df.shape[1] < n_expected_columns: - raise ValueError(f"The elements file has only {df.shape[1]} columns, but it should have at least {n_expected_columns}.") - - global_element_ids : NDArray[np.int_] = df.iloc[:, 0].to_numpy() - connectivity : NDArray[np.int_] = df.iloc[:, 3:6].to_numpy() + raise ValueError( + f"The elements file has only {df.shape[1]} columns, but it should have at least {n_expected_columns}." + ) + + global_element_ids: NDArray[np.int_] = df.iloc[:, 0].to_numpy() + connectivity: NDArray[np.int_] = df.iloc[:, 3:6].to_numpy() - return global_element_ids, connectivity \ No newline at end of file + return global_element_ids, connectivity diff --git a/src/reader.py b/src/reader.py index 4d2b7032..94dd59c2 100644 --- a/src/reader.py +++ b/src/reader.py @@ -32,29 +32,47 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: nodes_file: Path = mesh_root / f"part.{partition_id}.nodes" elements_file: Path = mesh_root / f"part.{partition_id}.elements" - assert header_file.is_file(), f"Header file {header_file} does not exist. Please ensure that this file exists in {mesh_root}." - assert nodes_file.is_file(), f"Nodes file {nodes_file} does not exist. Please ensure that this file exists in {mesh_root}." - assert elements_file.is_file(), f"Mesh file {elements_file} does not exist. Please ensure that this file exists in {mesh_root}." + assert ( + header_file.is_file() + ), f"Header file {header_file} does not exist. Please ensure that this file exists in {mesh_root}." + assert ( + nodes_file.is_file() + ), f"Nodes file {nodes_file} does not exist. Please ensure that this file exists in {mesh_root}." + assert ( + elements_file.is_file() + ), f"Mesh file {elements_file} does not exist. Please ensure that this file exists in {mesh_root}." # Parse header, nodes, and elements files n_vertices, n_cells = elmer.parser.parse_header(header_file) global_vertex_ids, x_vertices, y_vertices, z_vertices = elmer.parser.parse_nodes(nodes_file) local_vertex_ids = range(len(global_vertex_ids)) # use [0,1,...,n_vertices-1] to identify vertices locally - assert len(global_vertex_ids) == n_vertices, f"Number of vertex IDs in nodes file ({len(global_vertex_ids)}) does not match the header ({n_vertices})." - assert len(x_vertices) == n_vertices, f"Number of vertices in nodes file ({len(x_vertices)}) does not match the header ({n_vertices})." - assert len(y_vertices) == n_vertices, f"Number of vertices in nodes file ({len(y_vertices)}) does not match the header ({n_vertices})." - assert len(z_vertices) == n_vertices, f"Number of vertices in nodes file ({len(z_vertices)}) does not match the header ({n_vertices})." + assert ( + len(global_vertex_ids) == n_vertices + ), f"Number of vertex IDs in nodes file ({len(global_vertex_ids)}) does not match the header ({n_vertices})." + assert ( + len(x_vertices) == n_vertices + ), f"Number of vertices in nodes file ({len(x_vertices)}) does not match the header ({n_vertices})." + assert ( + len(y_vertices) == n_vertices + ), f"Number of vertices in nodes file ({len(y_vertices)}) does not match the header ({n_vertices})." + assert ( + len(z_vertices) == n_vertices + ), f"Number of vertices in nodes file ({len(z_vertices)}) does not match the header ({n_vertices})." global_cell_ids, global_cell_to_vertex = elmer.parser.parse_elements(elements_file) local_cell_ids = range(len(global_cell_ids)) # use [0,1,...,n_cells-1] to identify cells locally - vertex_l2g = {l:g for l,g in zip(local_vertex_ids, global_vertex_ids)} # map local ids [0,1,...n_verts] to global_vertex_ids - vertex_g2l = {g:l for l,g in vertex_l2g.items()} # invert dictionary to get map for global ids to local ids + vertex_l2g = { + l: g for l, g in zip(local_vertex_ids, global_vertex_ids) + } # map local ids [0,1,...n_verts] to global_vertex_ids + vertex_g2l = {g: l for l, g in vertex_l2g.items()} # invert dictionary to get map for global ids to local ids cell_to_vertex_local = np.array([[vertex_g2l[g_v] for g_v in c] for c in global_cell_to_vertex]) - assert len(global_cell_ids) == n_cells, f"Number of cell IDs in elements file ({len(global_cell_ids)}) does not match the header ({n_cells})." + assert ( + len(global_cell_ids) == n_cells + ), f"Number of cell IDs in elements file ({len(global_cell_ids)}) does not match the header ({n_cells})." return TriangleMesh( x_vertices=x_vertices, @@ -62,9 +80,10 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: z_vertices=z_vertices, cell_to_vertex=cell_to_vertex_local, vertex_ids=global_vertex_ids, - cell_ids=global_cell_ids + cell_ids=global_cell_ids, ) + def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]): """Read digital elevation model (DEM) file. Args: @@ -77,23 +96,29 @@ def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]): assert dem_file.is_file(), f"DEM file {dem_file} does not exist." import netCDF4 + nc = netCDF4.Dataset(dem_file) # find all matching ids in x direction - idx_x = np.searchsorted(nc['x'], xs) + idx_x = np.searchsorted(nc["x"], xs) # find all matching ids in y direction # nc['y'] is in reversed order. So finding the matches is a bit more complicated since np.searchsorted expects and array in ascending order # sorter=reversed(...) allows us to find indices in an array with descending order - idx_y = np.searchsorted(nc['y'],ys,sorter=list(reversed(range(len(nc['y']))))) - idx_y *= -1 # we have invert the sign since indices of the reversed array are counted from the back + idx_y = np.searchsorted(nc["y"], ys, sorter=list(reversed(range(len(nc["y"]))))) + idx_y *= -1 # we have invert the sign since indices of the reversed array are counted from the back idx_y += -1 - print(f"{np.sum(abs(xs - nc['x'][idx_x][:]) > 10e-5)} of {len(idx_x)} x-coordinates have a significant mismatch. surface value of a neighboring point will be used.") - print(f"{np.sum(abs(ys - nc['y'][idx_y][:]) > 10e-5)} of {len(idx_y)} y-coordinates have a significant mismatch. surface value of a neighboring point will be used.") + print( + f"{np.sum(abs(xs - nc['x'][idx_x][:]) > 10e-5)} of {len(idx_x)} x-coordinates have a significant mismatch. surface value of a neighboring point will be used." + ) + print( + f"{np.sum(abs(ys - nc['y'][idx_y][:]) > 10e-5)} of {len(idx_y)} y-coordinates have a significant mismatch. surface value of a neighboring point will be used." + ) - surf = nc['surface'][:] # get numpy array + surf = nc["surface"][:] # get numpy array return surf[idx_y, idx_x] # sample at given indexes and return + def read_matlab(mat_file: Path) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: """Read custom grid from a MATLAB file. Args: @@ -104,7 +129,9 @@ def read_matlab(mat_file: Path) -> tuple[NDArray[np.float64], NDArray[np.float64 raise Exception("Reading from MATLAB files is not implemented yet.") -def write_dem_as_elmer(mesh: TriangleMesh, h: NDArray[np.float64], dem_file: Path, allow_overwrite: bool = False) -> None: +def write_dem_as_elmer( + mesh: TriangleMesh, h: NDArray[np.float64], dem_file: Path, allow_overwrite: bool = False +) -> None: """Write digital elevation model to a file following the structure of an existing Elmer mesh Args: mesh (Mesh): The mesh object containing x and y vertices and vertex IDs. @@ -115,17 +142,14 @@ def write_dem_as_elmer(mesh: TriangleMesh, h: NDArray[np.float64], dem_file: Pat if not allow_overwrite: assert not dem_file.is_file(), f"DEM file {dem_file} already exists. Please choose a different file name." - assert len(h) == len(mesh.x_vertices), f"Height data length ({len(h)}) does not match number of vertices ({len(mesh.x_vertices)})." + assert len(h) == len( + mesh.x_vertices + ), f"Height data length ({len(h)}) does not match number of vertices ({len(mesh.x_vertices)})." import pandas as pd + # Create a DataFrame with the required structure - df = pd.DataFrame({ - 'Node ID': mesh.vertex_ids, - 'Node Type': -1, - 'x': mesh.x_vertices, - 'y': mesh.y_vertices, - 'z': h - }) + df = pd.DataFrame({"Node ID": mesh.vertex_ids, "Node Type": -1, "x": mesh.x_vertices, "y": mesh.y_vertices, "z": h}) def fortran_style_sci(x, precision=15): """Convert a number to Fortran-style scientific notation.""" @@ -133,32 +157,34 @@ def fortran_style_sci(x, precision=15): return f" 0.{''.join(['0']*precision)}E+00" # special case for zero; Fortran style requires leading space for positive numbers exp = int(np.floor(np.log10(abs(x)))) + 1 - mantissa = x / (10 ** exp) - sign = '-' if mantissa < 0 else ' ' + mantissa = x / (10**exp) + sign = "-" if mantissa < 0 else " " return sign + f"{abs(mantissa):.{precision}f}E{exp:+03d}" import csv + df.to_csv( dem_file, - sep=' ', + sep=" ", float_format=fortran_style_sci, index=False, header=False, quoting=csv.QUOTE_NONE, # do not use quotes around fields (avoids "..." around numbers if leading space is introduced for positive numbers) - escapechar='\\' # required when using QUOTE_NONE with special chars + escapechar="\\", # required when using QUOTE_NONE with special chars ) # Postprocess file to ensure it matches Elmer's expected format - with open(dem_file, 'r') as f: + with open(dem_file, "r") as f: content = ( - '\n'.join(line.rstrip() + ' ' for line in f.read().splitlines()) # append space to each line - .replace('\\ ', ' ') # replace escaped spaces with actual spaces - .replace(' -1 ', ' -1 ') # add another space after node type + "\n".join(line.rstrip() + " " for line in f.read().splitlines()) # append space to each line + .replace("\\ ", " ") # replace escaped spaces with actual spaces + .replace(" -1 ", " -1 ") # add another space after node type ) - with open(dem_file, 'w') as f: + with open(dem_file, "w") as f: f.write(content) + if __name__ == "__main__": parser = argparse.ArgumentParser(description="Read Elmer mesh and DEM files.") @@ -174,4 +200,4 @@ def fortran_style_sci(x, precision=15): y = mesh.y_vertices h = read_dem(args.dem, x, y) - write_dem_as_elmer(mesh, h, args.elmer_mesh / "MESH" / "mesh.nodes", allow_overwrite=True) \ No newline at end of file + write_dem_as_elmer(mesh, h, args.elmer_mesh / "MESH" / "mesh.nodes", allow_overwrite=True) diff --git a/src/utils.py b/src/utils.py index 032b9776..70717de6 100644 --- a/src/utils.py +++ b/src/utils.py @@ -3,22 +3,28 @@ import logging default_levels = { - 'file':logging.DEBUG, # log level for logging to file + "file": logging.DEBUG, # log level for logging to file 0: logging.INFO, # log level for rank 0 # 1: logging.DEBUG, # to log other ranks to console define log level here } + def setup_logging(log_levels=default_levels, file=None): root_logger = logging.getLogger() root_logger.setLevel(logging.DEBUG) - log_formatter = logging.Formatter(fmt=f'EBFM r{MPI.COMM_WORLD.rank} ' + '%(asctime)s.%(msecs)03d %(levelname)s %(module)s - %(funcName)s: %(message)s', - datefmt='%Y-%m-%d %H:%M:%S') + log_formatter = logging.Formatter( + fmt=f"EBFM r{MPI.COMM_WORLD.rank} " + + "%(asctime)s.%(msecs)03d %(levelname)s %(module)s - %(funcName)s: %(message)s", + datefmt="%Y-%m-%d %H:%M:%S", + ) if file: - file_handler = logging.FileHandler(filename=f'{file}' + f'.{MPI.COMM_WORLD.rank}' if MPI.COMM_WORLD.size > 1 else '') + file_handler = logging.FileHandler( + filename=f"{file}" + f".{MPI.COMM_WORLD.rank}" if MPI.COMM_WORLD.size > 1 else "" + ) file_handler.setFormatter(log_formatter) - file_handler.setLevel(log_levels['file']) + file_handler.setLevel(log_levels["file"]) root_logger.addHandler(file_handler) if MPI.COMM_WORLD.rank in log_levels.keys(): @@ -31,9 +37,8 @@ def setup_logging(log_levels=default_levels, file=None): stderr_handler = logging.StreamHandler(sys.stderr) stderr_handler.setFormatter(log_formatter) stderr_handler.setLevel(logging.ERROR) - root_logger.addHandler(stderr_handler) + root_logger.addHandler(stderr_handler) root_logger.info("Logging setup complete.") - return logging - + return logging From db12295a59f313ff53dcbcb91fe334b933e6439e Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 18:14:29 +0100 Subject: [PATCH 08/12] Manual (trivial) changes that flake8 flags. * black cannot fix this automatically. --- src/coupling.py | 56 +++++++++++++++++++++------ src/dummy.py | 37 +++++++++++++----- src/ebfm/FINAL_create_restart_file.py | 4 +- src/ebfm/INIT.py | 33 +++++++++++----- src/ebfm/LOOP_EBM.py | 9 ++++- src/ebfm/LOOP_SNOW.py | 6 +-- src/ebfm/LOOP_general_functions.py | 8 +++- src/ebfm/LOOP_write_to_file.py | 7 +++- src/elmer/parser.py | 14 ++++++- src/reader.py | 34 +++++++++++----- src/utils.py | 2 +- 11 files changed, 159 insertions(+), 51 deletions(-) diff --git a/src/coupling.py b/src/coupling.py index a779a5a7..89148496 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -85,7 +85,8 @@ class YACCoupler(Coupler): component: yac.Component = None grid: yac.UnstructuredGrid = None corner_points: yac.Points = None - # TODO use a single dictionary for both source and target fields; provide properties to create dicts containing source_fields and target_fields by comparing to exchange_type + # TODO use a single dictionary for both source and target fields; + # TODO provide properties to create dicts containing source_fields and target_fields by comparing to exchange_type source_fields: Dict[str, yac.Field] = {} target_fields: Dict[str, yac.Field] = {} @@ -154,9 +155,21 @@ def construct_coupling_elmer_ice(self, time: Dict[str, float]): name="T_ice", metadata="Near surface temperature at Ice surface (in K)", ), - FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="smb", metadata="??? (in ???)"), - FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="runoff", metadata="Runoff (in ???)"), - FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name="h", metadata="Surface height (in m)"), + FieldDefinition( + exchange_type=yac.ExchangeType.SOURCE, + name="smb", + metadata="??? (in ???)", + ), + FieldDefinition( + exchange_type=yac.ExchangeType.SOURCE, + name="runoff", + metadata="Runoff (in ???)", + ), + FieldDefinition( + exchange_type=yac.ExchangeType.TARGET, + name="h", + metadata="Surface height (in m)", + ), ] timestep_value = days_to_iso(time["dt"]) @@ -175,7 +188,10 @@ def construct_coupling_elmer_ice(self, time: Dict[str, float]): timestep.format, ) self.interface.def_field_metadata( - field.component_name, field.grid_name, field_definition.name, field_definition.metadata.encode("utf-8") + field.component_name, + field.grid_name, + field_definition.name, + field_definition.metadata.encode("utf-8"), ) if field_definition.exchange_type == yac.ExchangeType.SOURCE: @@ -190,9 +206,15 @@ def construct_coupling_icon_atmo(self, time: Dict[str, float]): FieldDefinition = namedtuple("FieldDefinition", ["exchange_type", "name", "metadata"]) # TODO: Get hard-coded data below from dummies/EBFM/ebfm-config.yaml field_definitions = [ - # FieldDefinition(exchange_type=yac.ExchangeType.SOURCE, name="albedo", metadata="Albedo of the ice surface (in ???)"), + # FieldDefinition( + # exchange_type=yac.ExchangeType.SOURCE, + # name="albedo", + # metadata="Albedo of the ice surface (in ???)" + # ), FieldDefinition( - exchange_type=yac.ExchangeType.TARGET, name="pr", metadata="Precipitation rate (in kg m-2 s-1)" + exchange_type=yac.ExchangeType.TARGET, + name="pr", + metadata="Precipitation rate (in kg m-2 s-1)", ), FieldDefinition( exchange_type=yac.ExchangeType.TARGET, @@ -216,8 +238,16 @@ def construct_coupling_icon_atmo(self, time: Dict[str, float]): FieldDefinition( exchange_type=yac.ExchangeType.TARGET, name="tas", metadata="Temperature at surface (in K)" ), - # FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='huss', metadata='Specific humidity at surface (in kg kg-1)'), - # FieldDefinition(exchange_type=yac.ExchangeType.TARGET, name='sfcPressure', metadata='Surface pressure (in Pa)') + # FieldDefinition( + # exchange_type=yac.ExchangeType.TARGET, + # name="huss", + # metadata="Specific humidity at surface (in kg kg-1)" + # ), + # FieldDefinition( + # exchange_type=yac.ExchangeType.TARGET, + # name="sfcPressure", + # metadata="Surface pressure (in Pa)" + # ), ] timestep_value = days_to_iso(time["dt"]) @@ -250,9 +280,11 @@ def construct_coupling_post_sync(self): for field in self.target_fields.values(): is_defined = self.interface.get_field_is_defined(field.component_name, field.grid_name, field.name) - assert ( - is_defined - ), f"Field '{field.name}' is not defined in YAC for component '{field.component_name}' and grid '{field.grid_name}'." + assert is_defined, ( + f"Field '{field.name}' is not defined in YAC for component '{field.component_name}' and " + f"grid '{field.grid_name}'." + ) + field_role = self.interface.get_field_role(field.component_name, field.grid_name, field.name) assert ( field_role == yac.ExchangeType.TARGET diff --git a/src/dummy.py b/src/dummy.py index 8e635ed1..49f2cf4d 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -2,7 +2,14 @@ from pathlib import Path import argparse -from ebfm import INIT, LOOP_general_functions, LOOP_climate_forcing, LOOP_EBM, LOOP_SNOW, LOOP_mass_balance +from ebfm import ( + INIT, + LOOP_general_functions, + LOOP_climate_forcing, + LOOP_EBM, + LOOP_SNOW, + LOOP_mass_balance, +) from ebfm import LOOP_write_to_file, FINAL_create_restart_file from ebfm.grid import GridInputType @@ -30,17 +37,27 @@ def main(): parser = argparse.ArgumentParser() parser.add_argument( - "--couple-to-elmer-ice", action="store_true", help="Enable coupling with Elmer/Ice models via YAC" + "--couple-to-elmer-ice", + action="store_true", + help="Enable coupling with Elmer/Ice models via YAC", ) - parser.add_argument("--couple-to-icon-atmo", action="store_true", help="Enable coupling with ICON via YAC") + parser.add_argument( + "--couple-to-icon-atmo", + action="store_true", + help="Enable coupling with ICON via YAC", + ) parser.add_argument( - "--coupler-config", type=Path, help="Path to the coupling configuration file (YAC coupler_config.yaml)." + "--coupler-config", + type=Path, + help="Path to the coupling configuration file (YAC coupler_config.yaml).", ) parser.add_argument( - "--elmer-mesh", type=Path, help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required." + "--elmer-mesh", + type=Path, + help="Path to the Elmer mesh file. Either --elmer-mesh or --matlab-mesh is required.", ) parser.add_argument( @@ -52,7 +69,8 @@ def main(): parser.add_argument( "--netcdf-mesh", type=Path, - help="Path to the NetCDF mesh file. Optional if using --elmer-mesh. If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.", + help="Path to the NetCDF mesh file. Optional if using --elmer-mesh." + " If --netcdf-mesh is provided elevations will be read from the given NetCDF mesh file.", ) parser.add_argument( @@ -65,13 +83,14 @@ def main(): "--use-part", type=int, default=MPI.COMM_WORLD.rank + 1, - help="If using a partitioned Elmer mesh, allows to specify which partition ID to use for this run. If not provided, the MPI rank + 1 will be used as partition ID.", + help="If using a partitioned Elmer mesh, allows to specify which partition ID to use for this run. " + "If not provided, the MPI rank + 1 will be used as partition ID.", ) args = parser.parse_args() - logger.info(f"Done parsing command line arguments.") - logger.debug(f"Parsed the following command line arguments:") + logger.info("Done parsing command line arguments.") + logger.debug("Parsed the following command line arguments:") for arg, val in vars(args).items(): logger.debug(f" {arg}: {val}") diff --git a/src/ebfm/FINAL_create_restart_file.py b/src/ebfm/FINAL_create_restart_file.py index c571b1f9..27e147de 100644 --- a/src/ebfm/FINAL_create_restart_file.py +++ b/src/ebfm/FINAL_create_restart_file.py @@ -53,7 +53,9 @@ def create_boot_file(): return True OUT["timelastsnow_netCDF"] = date2num( - OUT["timelastsnow"], units="days since 1970-01-01 00:00:00", calendar="gregorian" + OUT["timelastsnow"], + units="days since 1970-01-01 00:00:00", + calendar="gregorian", ) OUT = { "subZ": OUT["subZ"], diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index aa7941e6..8e4bae33 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -172,9 +172,10 @@ def init_grid(grid, io, args: Namespace): grid["is_partitioned"] = args.is_partitioned_elmer_mesh if grid["is_partitioned"]: - assert ( - args.netcdf_mesh - ), "--is-partitioned-elmer-mesh requires --netcdf-mesh. (Without --netcdf-mesh should also work but is untested.)" + assert args.netcdf_mesh, ( + "--is-partitioned-elmer-mesh requires --netcdf-mesh. " + "(Without --netcdf-mesh should also work but is untested.)" + ) logger.info("Using partitioned grid...") else: logger.info("Using non-partitioned grid...") @@ -187,13 +188,18 @@ def init_grid(grid, io, args: Namespace): grid["input_type"] = GridInputType.ELMER else: logger.error( - f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. Please please refer to the documentation for correct configuration." + f"Invalid grid configuration. EBFM supports the grid types {[t.name for t in GridInputType]}. " + "Please please refer to the documentation for correct configuration." ) raise Exception("Invalid grid configuration.") if grid["input_type"] is GridInputType.CUSTOM: # Read grid from Elmer, elevations from BedMachine if grid["is_partitioned"]: - mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh, is_partitioned=True, partition_id=args.use_part) + mesh: Mesh = read_elmer_mesh( + mesh_root=args.elmer_mesh, + is_partitioned=True, + partition_id=args.use_part, + ) else: mesh: Mesh = read_elmer_mesh(mesh_root=args.elmer_mesh) @@ -211,10 +217,15 @@ def init_grid(grid, io, args: Namespace): # TODO later add slope # dzdx, dzdy = mesh.dzdy, mesh.dzdy elif grid["input_type"] is GridInputType.ELMER: # Read grid and elevations from Elmer - mesh: Mesh = read_elmer_mesh( - args.elmer_mesh - ) # assuming mesh/MESH/mesh.nodes contains DEM data in the z component; see mesh/README.md for the required preprocessing steps. - grid["x"], grid["y"], grid["z"] = mesh.x_vertices, mesh.y_vertices, mesh.z_vertices + mesh: Mesh = read_elmer_mesh(args.elmer_mesh) + + # assuming mesh/MESH/mesh.nodes contains DEM data in the z component + # see mesh/README.md for the required preprocessing steps. + grid["x"], grid["y"], grid["z"] = ( + mesh.x_vertices, + mesh.y_vertices, + mesh.z_vertices, + ) grid["z"] = np.random.uniform(0, 100, size=len(grid["x"])) # test values! grid["slope_x"] = np.zeros_like(grid["x"]) # test values! grid["slope_y"] = np.zeros_like(grid["x"]) # test values! @@ -414,7 +425,9 @@ def init_initial_conditions(C, grid, io, time): OUT[var_name] = var_data OUT["timelastsnow"] = num2date( - OUT["timelastsnow_netCDF"], units="days since 1970-01-01 00:00:00", calendar="gregorian" + OUT["timelastsnow_netCDF"], + units="days since 1970-01-01 00:00:00", + calendar="gregorian", ) else: diff --git a/src/ebfm/LOOP_EBM.py b/src/ebfm/LOOP_EBM.py index 4937afab..8d4c6dd9 100644 --- a/src/ebfm/LOOP_EBM.py +++ b/src/ebfm/LOOP_EBM.py @@ -1,6 +1,13 @@ import numpy as np -from ebfm import LOOP_EBM_SHF, LOOP_EBM_GHF, LOOP_EBM_LHF, LOOP_EBM_LWin, LOOP_EBM_LWout, LOOP_EBM_SWin +from ebfm import ( + LOOP_EBM_SHF, + LOOP_EBM_GHF, + LOOP_EBM_LHF, + LOOP_EBM_LWin, + LOOP_EBM_LWout, + LOOP_EBM_SWin, +) from ebfm import LOOP_EBM_SWout, LOOP_EBM_insolation from coupling import Coupler diff --git a/src/ebfm/LOOP_SNOW.py b/src/ebfm/LOOP_SNOW.py index b95b0efc..fe7e65a7 100644 --- a/src/ebfm/LOOP_SNOW.py +++ b/src/ebfm/LOOP_SNOW.py @@ -560,9 +560,9 @@ def percolation_refreezing_and_storage(): # Update density after refreezing OUT["subD"] += RS / OUT["subZ"] - ############################################################# - ### REFREEZING OF IRREDUCIBLE WATER - ############################################################# + ########################################################### + # REFREEZING OF IRREDUCIBLE WATER + ########################################################### # Determine whether cold content or density limits the amount of refreezing OUT["cpi"] = 152.2 + 7.122 * OUT["subT"] c1 = OUT["cpi"] * OUT["subD"] * OUT["subZ"] * (C["T0"] - OUT["subT"]) / C["Lm"] diff --git a/src/ebfm/LOOP_general_functions.py b/src/ebfm/LOOP_general_functions.py index e4dec768..377db61e 100644 --- a/src/ebfm/LOOP_general_functions.py +++ b/src/ebfm/LOOP_general_functions.py @@ -9,7 +9,8 @@ def print_time(t, time): Parameters: t (int): Current time step. - time (dict): Dictionary containing time-related parameters, such as start time (ts), end time (te), and time step (dt). + time (dict): Dictionary containing time-related parameters, such as start time (ts), end time (te), and time + step (dt). Returns: dict: Updated time dictionary with calculated current and previous times. @@ -27,7 +28,10 @@ def print_time(t, time): # Print the current time parameters print( - f"Year: {tempdate.tm_year:4d} Month: {tempdate.tm_mon:2d} Day: {tempdate.tm_mday:2d} Hour: {tempdate.tm_hour:2d}" + f"Year: {tempdate.tm_year:4d} " + f"Month: {tempdate.tm_mon:2d} " + f"Day: {tempdate.tm_mday:2d} " + f"Hour: {tempdate.tm_hour:2d}" ) return time diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index a30ca699..ff5d71d6 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -162,7 +162,12 @@ def save_netCDF_file(): zlib=True, complevel=4, fill_value=-9999.0, # Fill missing values - chunksizes=(1, grid["x_2D"].shape[0], grid["x_2D"].shape[1], grid["nl"]), + chunksizes=( + 1, + grid["x_2D"].shape[0], + grid["x_2D"].shape[1], + grid["nl"], + ), ) else: # Define variable as 3D: (time, y, x) diff --git a/src/elmer/parser.py b/src/elmer/parser.py index 163e835d..bb88ade3 100644 --- a/src/elmer/parser.py +++ b/src/elmer/parser.py @@ -55,7 +55,12 @@ def parse_nodes( Returns: tuple: Global node IDs, x coordinates, y coordinates, and z coordinates. """ - df = pd.read_csv(nodes_file, sep="\\s+", header=None, dtype={0: int, 2: float, 3: float, 4: float}) + df = pd.read_csv( + nodes_file, + sep="\\s+", + header=None, + dtype={0: int, 2: float, 3: float, 4: float}, + ) # A valid mesh.nodes file in Elmer should have exactly 5 columns: # 1. Global Node ID @@ -84,7 +89,12 @@ def parse_elements(elements_file: Path) -> tuple[NDArray[np.int_], NDArray[np.in Returns: tuple: Element IDs and connectivity. """ - df = pd.read_csv(elements_file, sep="\\s+", header=None, dtype={0: int, 3: float, 4: float, 5: float}) + df = pd.read_csv( + elements_file, + sep="\\s+", + header=None, + dtype={0: int, 3: float, 4: float, 5: float}, + ) # A valid mesh.elements file in Elmer should have at least 2 columns: # 1. Global Element ID diff --git a/src/reader.py b/src/reader.py index 94dd59c2..d21a7307 100644 --- a/src/reader.py +++ b/src/reader.py @@ -102,24 +102,28 @@ def read_dem(dem_file: Path, xs: NDArray[np.float64], ys: NDArray[np.float64]): idx_x = np.searchsorted(nc["x"], xs) # find all matching ids in y direction - # nc['y'] is in reversed order. So finding the matches is a bit more complicated since np.searchsorted expects and array in ascending order - # sorter=reversed(...) allows us to find indices in an array with descending order + # nc['y'] is in reversed order. So finding the matches is a bit more complicated since np.searchsorted expects and + # array in ascending order sorter=reversed(...) allows us to find indices in an array with descending order idx_y = np.searchsorted(nc["y"], ys, sorter=list(reversed(range(len(nc["y"]))))) idx_y *= -1 # we have invert the sign since indices of the reversed array are counted from the back idx_y += -1 print( - f"{np.sum(abs(xs - nc['x'][idx_x][:]) > 10e-5)} of {len(idx_x)} x-coordinates have a significant mismatch. surface value of a neighboring point will be used." + f"{np.sum(abs(xs - nc['x'][idx_x][:]) > 10e-5)} of {len(idx_x)} x-coordinates have a significant mismatch." + " Surface value of a neighboring point will be used." ) print( - f"{np.sum(abs(ys - nc['y'][idx_y][:]) > 10e-5)} of {len(idx_y)} y-coordinates have a significant mismatch. surface value of a neighboring point will be used." + f"{np.sum(abs(ys - nc['y'][idx_y][:]) > 10e-5)} of {len(idx_y)} y-coordinates have a significant mismatch." + " Surface value of a neighboring point will be used." ) surf = nc["surface"][:] # get numpy array return surf[idx_y, idx_x] # sample at given indexes and return -def read_matlab(mat_file: Path) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: +def read_matlab( + mat_file: Path, +) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: """Read custom grid from a MATLAB file. Args: mat_file (Path): Path to the MATLAB file. @@ -130,7 +134,10 @@ def read_matlab(mat_file: Path) -> tuple[NDArray[np.float64], NDArray[np.float64 def write_dem_as_elmer( - mesh: TriangleMesh, h: NDArray[np.float64], dem_file: Path, allow_overwrite: bool = False + mesh: TriangleMesh, + h: NDArray[np.float64], + dem_file: Path, + allow_overwrite: bool = False, ) -> None: """Write digital elevation model to a file following the structure of an existing Elmer mesh Args: @@ -149,12 +156,21 @@ def write_dem_as_elmer( import pandas as pd # Create a DataFrame with the required structure - df = pd.DataFrame({"Node ID": mesh.vertex_ids, "Node Type": -1, "x": mesh.x_vertices, "y": mesh.y_vertices, "z": h}) + df = pd.DataFrame( + { + "Node ID": mesh.vertex_ids, + "Node Type": -1, + "x": mesh.x_vertices, + "y": mesh.y_vertices, + "z": h, + } + ) def fortran_style_sci(x, precision=15): """Convert a number to Fortran-style scientific notation.""" if x == 0: - return f" 0.{''.join(['0']*precision)}E+00" # special case for zero; Fortran style requires leading space for positive numbers + # special case for zero; Fortran style requires leading space for positive numbers + return f" 0.{''.join(['0'] * precision)}E+00" exp = int(np.floor(np.log10(abs(x)))) + 1 mantissa = x / (10**exp) @@ -169,8 +185,8 @@ def fortran_style_sci(x, precision=15): float_format=fortran_style_sci, index=False, header=False, - quoting=csv.QUOTE_NONE, # do not use quotes around fields (avoids "..." around numbers if leading space is introduced for positive numbers) escapechar="\\", # required when using QUOTE_NONE with special chars + quoting=csv.QUOTE_NONE, # do not use quotes around fields ) # Postprocess file to ensure it matches Elmer's expected format diff --git a/src/utils.py b/src/utils.py index 70717de6..10354d52 100644 --- a/src/utils.py +++ b/src/utils.py @@ -21,7 +21,7 @@ def setup_logging(log_levels=default_levels, file=None): if file: file_handler = logging.FileHandler( - filename=f"{file}" + f".{MPI.COMM_WORLD.rank}" if MPI.COMM_WORLD.size > 1 else "" + filename=(f"{file}" + f".{MPI.COMM_WORLD.rank}" if MPI.COMM_WORLD.size > 1 else "") ) file_handler.setFormatter(log_formatter) file_handler.setLevel(log_levels["file"]) From 3eaa03619943628ef2096200cd8e8484f2b38cda Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 18:15:01 +0100 Subject: [PATCH 09/12] Ignore some rules. --- .flake8 | 1 + 1 file changed, 1 insertion(+) diff --git a/.flake8 b/.flake8 index aa079ec5..52c61bc8 100644 --- a/.flake8 +++ b/.flake8 @@ -1,2 +1,3 @@ [flake8] +ignore = E203, W503 max-line-length=120 From 239d3d1ccf45154b6c80d3d67af950382e40d029 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 18:16:10 +0100 Subject: [PATCH 10/12] Remove unused imports. --- src/coupling.py | 3 +-- src/dummy.py | 3 --- src/ebfm/INIT.py | 2 -- src/ebfm/LOOP_EBM_GHF.py | 3 --- src/ebfm/LOOP_EBM_LWin.py | 3 --- src/ebfm/LOOP_EBM_LWout.py | 3 --- src/ebfm/LOOP_general_functions.py | 4 +--- src/ebfm/LOOP_write_to_file.py | 1 - 8 files changed, 2 insertions(+), 20 deletions(-) diff --git a/src/coupling.py b/src/coupling.py index 89148496..7907807c 100644 --- a/src/coupling.py +++ b/src/coupling.py @@ -1,9 +1,8 @@ import yac from pathlib import Path -import pyproj import numpy as np from collections import namedtuple -from typing import List, Dict +from typing import Dict from elmer.mesh import Mesh as Grid # for now use an alias diff --git a/src/dummy.py b/src/dummy.py index 49f2cf4d..d9263889 100644 --- a/src/dummy.py +++ b/src/dummy.py @@ -17,9 +17,6 @@ from utils import setup_logging import logging -from typing import NamedTuple - - log_levels = { "file": logging.DEBUG, # log level for logging to file 0: logging.INFO, # log level for rank 0 diff --git a/src/ebfm/INIT.py b/src/ebfm/INIT.py index 8e4bae33..1917c190 100644 --- a/src/ebfm/INIT.py +++ b/src/ebfm/INIT.py @@ -1,6 +1,5 @@ import os from datetime import datetime -from enum import Enum from argparse import Namespace from typing import Any import numpy as np @@ -14,7 +13,6 @@ from elmer.mesh import Mesh from ebfm.grid import GridInputType -from mpi4py import MPI import logging diff --git a/src/ebfm/LOOP_EBM_GHF.py b/src/ebfm/LOOP_EBM_GHF.py index ccccf72c..47d9b6c6 100644 --- a/src/ebfm/LOOP_EBM_GHF.py +++ b/src/ebfm/LOOP_EBM_GHF.py @@ -1,6 +1,3 @@ -import numpy as np - - def main(Tsurf, OUT, cond, GHF_k, GHF_C): """ Calculates the subsurface heat flux (GHF) based on effective conductivity diff --git a/src/ebfm/LOOP_EBM_LWin.py b/src/ebfm/LOOP_EBM_LWin.py index 17f82caa..b0ddb2db 100644 --- a/src/ebfm/LOOP_EBM_LWin.py +++ b/src/ebfm/LOOP_EBM_LWin.py @@ -1,6 +1,3 @@ -import numpy as np - - def main(C, IN): """ Computes incoming longwave radiation (LWin) based on atmospheric emissivity. diff --git a/src/ebfm/LOOP_EBM_LWout.py b/src/ebfm/LOOP_EBM_LWout.py index d97b3ad0..402b9895 100644 --- a/src/ebfm/LOOP_EBM_LWout.py +++ b/src/ebfm/LOOP_EBM_LWout.py @@ -1,6 +1,3 @@ -import numpy as np - - def main(C, Tsurf): """ Calculates the outgoing longwave radiation (LWout) based on blackbody diff --git a/src/ebfm/LOOP_general_functions.py b/src/ebfm/LOOP_general_functions.py index 377db61e..176c7e7f 100644 --- a/src/ebfm/LOOP_general_functions.py +++ b/src/ebfm/LOOP_general_functions.py @@ -1,6 +1,4 @@ -from datetime import datetime, timedelta - -from netCDF4 import date2num +from datetime import timedelta def print_time(t, time): diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index ff5d71d6..037246cf 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -1,4 +1,3 @@ -import datetime import os import numpy as np from netCDF4 import Dataset, date2num From 350c7e54d42ebb576779463fe8c3ec56b2442fb4 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 18:17:55 +0100 Subject: [PATCH 11/12] Remove dead code and unused variables. --- src/ebfm/LOOP_EBM_insolation.py | 3 --- src/ebfm/LOOP_SNOW.py | 1 - src/ebfm/LOOP_general_functions.py | 4 ---- src/ebfm/LOOP_write_to_file.py | 5 ++--- src/reader.py | 1 - 5 files changed, 2 insertions(+), 12 deletions(-) diff --git a/src/ebfm/LOOP_EBM_insolation.py b/src/ebfm/LOOP_EBM_insolation.py index 4243ed15..d0a4c19c 100644 --- a/src/ebfm/LOOP_EBM_insolation.py +++ b/src/ebfm/LOOP_EBM_insolation.py @@ -105,9 +105,6 @@ def main(grid, time2, OUT): ) ddy[azimuth > 0.75 * np.pi] = 1 - # Initialize the shade_2D array - mask = grid["mask_2D"] - # Shading routine z_flat = grid["z_2D"].ravel() z_current = z_flat[grid["ind"]] diff --git a/src/ebfm/LOOP_SNOW.py b/src/ebfm/LOOP_SNOW.py index fe7e65a7..8f819f27 100644 --- a/src/ebfm/LOOP_SNOW.py +++ b/src/ebfm/LOOP_SNOW.py @@ -331,7 +331,6 @@ def heat_conduction(): """ Calculate heat diffusion and update temperatures """ - nl = grid["nl"] dz1 = (OUT["subZ"][:, 0] + 0.5 * OUT["subZ"][:, 1]) ** 2 dz2 = 0.5 * (OUT["subZ"][:, 2:] + OUT["subZ"][:, 1:-1]) ** 2 kk = 0.138 - 1.01e-3 * OUT["subD"] + 3.233e-6 * OUT["subD"] ** 2 # Effective conductivity diff --git a/src/ebfm/LOOP_general_functions.py b/src/ebfm/LOOP_general_functions.py index 176c7e7f..5cdc7998 100644 --- a/src/ebfm/LOOP_general_functions.py +++ b/src/ebfm/LOOP_general_functions.py @@ -13,10 +13,6 @@ def print_time(t, time): Returns: dict: Updated time dictionary with calculated current and previous times. """ - # Convert time parameters to datetime objects - # time['TS'] = datetime.strptime(time['ts'], '%Y-%m-%d %H:%M:%S') # Start time - # time['TE'] = datetime.strptime(time['te'], '%Y-%m-%d %H:%M:%S') # End time - # Calculate current and previous model times based on the time step time["TCUR"] = time["ts"] + timedelta(days=(t - 1) * time["dt"]) time["TPREV"] = time["ts"] + timedelta(days=(t - 2) * time["dt"]) diff --git a/src/ebfm/LOOP_write_to_file.py b/src/ebfm/LOOP_write_to_file.py index 037246cf..24eead72 100644 --- a/src/ebfm/LOOP_write_to_file.py +++ b/src/ebfm/LOOP_write_to_file.py @@ -104,10 +104,9 @@ def save_binary_files(): for file in io["fid"].values(): file.close() + # TODO: WORK IN PROGRESS # Prepare the runinfo dictionary - runinfo = {"grid": grid, "time": time, "IOout": io, "Cout": C} - - ################ WORK IN PROGRESS ############### + # runinfo = {"grid": grid, "time": time, "IOout": io, "Cout": C} return True diff --git a/src/reader.py b/src/reader.py index d21a7307..3eb8f56b 100644 --- a/src/reader.py +++ b/src/reader.py @@ -61,7 +61,6 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: ), f"Number of vertices in nodes file ({len(z_vertices)}) does not match the header ({n_vertices})." global_cell_ids, global_cell_to_vertex = elmer.parser.parse_elements(elements_file) - local_cell_ids = range(len(global_cell_ids)) # use [0,1,...,n_cells-1] to identify cells locally vertex_l2g = { l: g for l, g in zip(local_vertex_ids, global_vertex_ids) From ab8559260e2fb453ba3e30b7aecdd72992ce6d07 Mon Sep 17 00:00:00 2001 From: Benjamin Rodenberg Date: Wed, 5 Nov 2025 18:18:33 +0100 Subject: [PATCH 12/12] Replace l with loc, See flake8 rule E741. --- src/reader.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/reader.py b/src/reader.py index 3eb8f56b..d5bf2cc7 100644 --- a/src/reader.py +++ b/src/reader.py @@ -63,9 +63,11 @@ def read_elmer_mesh(mesh_root: Path, is_partitioned: bool = False, partition_id: global_cell_ids, global_cell_to_vertex = elmer.parser.parse_elements(elements_file) vertex_l2g = { - l: g for l, g in zip(local_vertex_ids, global_vertex_ids) + loc: glob for loc, glob in zip(local_vertex_ids, global_vertex_ids) } # map local ids [0,1,...n_verts] to global_vertex_ids - vertex_g2l = {g: l for l, g in vertex_l2g.items()} # invert dictionary to get map for global ids to local ids + vertex_g2l = { + glob: loc for loc, glob in vertex_l2g.items() + } # invert dictionary to get map for global ids to local ids cell_to_vertex_local = np.array([[vertex_g2l[g_v] for g_v in c] for c in global_cell_to_vertex])