diff --git a/conda_requirements.txt b/conda_requirements.txt
index a1caeb7..c2ce2ed 100644
--- a/conda_requirements.txt
+++ b/conda_requirements.txt
@@ -15,6 +15,7 @@ nco
 cdo
 cartopy
 vtk
+pydantic
 richdem
 coveralls
 codecov
diff --git a/cosipy/config.py b/cosipy/config.py
index d8944d2..1108986 100644
--- a/cosipy/config.py
+++ b/cosipy/config.py
@@ -6,6 +6,12 @@
 import sys
 from importlib.metadata import entry_points
 from pathlib import Path
+from typing import Annotated, Literal, Optional
+
+from pydantic import BaseModel, ConfigDict, Field, model_validator
+from pydantic.types import StringConstraints
+from rich import print  # noqa: A004
+from typing_extensions import Self
 
 if sys.version_info >= (3, 11):
     import tomllib
@@ -95,6 +101,74 @@ def get_user_arguments() -> argparse.Namespace:
     return arguments
 
 
+DatetimeStr = Annotated[
+    str,
+    StringConstraints(
+        strip_whitespace=True, pattern=r"\d{4}-[01]\d-[0-3]\dT[0-2]\d:[0-5]\d"
+    ),
+]
+
+
+class CosipyConfigModel(BaseModel):
+    """COSIPY configuration model."""
+
+    model_config = ConfigDict(from_attributes=True)
+    time_start: DatetimeStr = Field(
+        description="Start time of the simulation in ISO format"
+    )
+    time_end: DatetimeStr = Field(description="End time of the simulation in ISO format")
+    data_path: Path = Field(description="Path to the data directory")
+    input_netcdf: Path = Field(description="Input NetCDF file path")
+    output_prefix: str = Field(description="Prefix for output files")
+    restart: bool = Field(description="Restart flag")
+    stake_evaluation: bool = Field(description="Flag for stake data evaluation")
+    stakes_loc_file: Path = Field(description="Path to stake location file")
+    stakes_data_file: Path = Field(description="Path to stake data file")
+    eval_method: Literal["rmse"] = Field(
+        "rmse", description="Evaluation method for simulations"
+    )
+    obs_type: Literal["mb", "snowheight"] = Field(description="Type of stake data used")
+    WRF: bool = Field(description="Flag for WRF input")
+    WRF_X_CSPY: bool = Field(description="Interactive simulation with WRF flag")
+    northing: str = Field(description="Name of northing dimension")
+    easting: str = Field(description="Name of easting dimension")
+    compression_level: int = Field(
+        ge=0, le=9, description="Output NetCDF compression level"
+    )
+    slurm_use: bool = Field(description="Use SLURM flag")
+    workers: Optional[int] = Field(
+        default=None,
+        ge=0,
+        description="""
+        Setting is only used is slurm_use is False.
+        Number of workers (cores), with 0 all available cores are used.
+        """,
+    )
+    local_port: int = Field(default=8786, gt=0, description="Port for local cluster")
+    full_field: bool = Field(description="Flag for writing full fields to file")
+    force_use_TP: bool = Field(..., description="Total precipitation flag")
+    force_use_N: bool = Field(..., description="Cloud cover fraction flag")
+    tile: bool = Field(description="Flag for tiling")
+    xstart: int = Field(ge=0, description="Start x index")
+    xend: int = Field(ge=0, description="End x index")
+    ystart: int = Field(ge=0, description="Start y index")
+    yend: int = Field(ge=0, description="End y index")
+    output_atm: str = Field(description="Atmospheric output variables")
+    output_internal: str = Field(description="Internal output variables")
+    output_full: str = Field(description="Full output variables")
+
+    @model_validator(mode="after")
+    def validate_output_variables(self) -> Self:
+        if self.WRF:
+            self.northing = "south_north"
+            self.easting = "west_east"
+        if self.WRF_X_CSPY:
+            self.full_field = True
+        if self.workers == 0:
+            self.workers = None
+        return self
+
+
 def get_help():
     """Print help for commands."""
     parser = set_parser()
@@ -135,7 +209,7 @@ class TomlLoader(object):
     """Load and parse configuration files."""
 
     @staticmethod
-    def get_raw_toml(file_path: str = "./config.toml") -> dict:
+    def get_raw_toml(file_path: Path = default_path) -> dict:
         """Open and load .toml configuration file.
 
         Args:
@@ -144,21 +218,20 @@ def get_raw_toml(file_path: str = "./config.toml") -> dict:
         Returns:
             Loaded .toml data.
         """
-        with open(file_path, "rb") as f:
-            raw_config = tomllib.load(f)
-
-        return raw_config
+        with file_path.open("rb") as f:
+            return tomllib.load(f)
 
-    @classmethod
-    def set_config_values(cls, config_table: dict):
+    @staticmethod
+    def flatten(config_table: dict[str, dict]) -> dict:
         """Overwrite attributes with configuration data.
 
         Args:
             config_table: Loaded .toml data.
         """
-        for _, table in config_table.items():
-            for k, v in table.items():
-                setattr(cls, k, v)
+        flat_dict = {}
+        for table in config_table.values():
+            flat_dict = {**flat_dict, **table}
+        return flat_dict
 
 
 class Config(TomlLoader):
@@ -168,37 +241,45 @@ class Config(TomlLoader):
     .toml file.
     """
 
-    def __init__(self):
-        self.args = get_user_arguments()
-        self.load(self.args.config_path)
-
-    @classmethod
-    def load(cls, path: str = "./config.toml"):
-        raw_toml = cls.get_raw_toml(path)
-        parsed_toml = cls.set_correct_config(raw_toml)
-        cls.set_config_values(parsed_toml)
-
-    @classmethod
-    def set_correct_config(cls, config_table: dict) -> dict:
-        """Adjust invalid or mutually exclusive configuration values.
+    def __init__(self, path: Path = default_path) -> None:
+        raw_toml = self.get_raw_toml(path)
+        self.raw_toml = self.flatten(raw_toml)
 
-        Args:
-            config_table: Loaded .toml data.
+    def validate(self) -> CosipyConfigModel:
+        """Validate configuration using Pydantic class.
 
         Returns:
-            Adjusted .toml data.
+            CosipyConfigModel: Validated configuration.
         """
-        # WRF Compatibility
-        if config_table["DIMENSIONS"]["WRF"]:
-            config_table["DIMENSIONS"]["northing"] = "south_north"
-            config_table["DIMENSIONS"]["easting"] = "west_east"
-        if config_table["DIMENSIONS"]["WRF_X_CSPY"]:
-            config_table["FULL_FIELDS"]["full_field"] = True
-        # TOML doesn't support null values
-        if config_table["PARALLELIZATION"]["workers"] == 0:
-            config_table["PARALLELIZATION"]["workers"] = None
+        return CosipyConfigModel(**self.raw_toml)
+
+
+ShebangStr = Annotated[str, StringConstraints(strip_whitespace=True, pattern=r"^#!")]
+
+
+class SlurmConfigModel(BaseModel):
+    """Slurm configuration model."""
+
+    account: str = Field(description="Slurm account/group")
+    name: str = Field(description="Equivalent to Slurm parameter `--job-name`")
+    queue: str = Field(description="Queue name")
+    slurm_parameters: list[str] = Field(description="Additional Slurm parameters")
+    shebang: ShebangStr = Field(description="Shebang string")
+    local_directory: Path = Field(description="Local directory")
+    port: int = Field(description="Network port number")
+    cores: int = Field(description="One grid point per core")
+    nodes: int = Field(description="Grid points submitted in one sbatch script")
+    processes: int = Field(description="Number of processes")
+    memory: str = Field(description="Memory per process")
+    memory_per_process: Optional[int] = Field(gt=0, description="Memory per process")
 
-        return config_table
+    @model_validator(mode="after")
+    def validate_output_variables(self):
+        if self.memory_per_process:
+            memory = self.memory_per_process * self.cores
+            self.memory = f"{memory}GB"
+
+        return self
 
 
 class SlurmConfig(TomlLoader):
@@ -222,43 +303,27 @@ class SlurmConfig(TomlLoader):
         slurm_parameters (List[str]): Additional Slurm parameters.
     """
 
-    def __init__(self):
-        self.args = get_user_arguments()
-        self.load(self.args.slurm_path)
-
-    @classmethod
-    def load(cls, path: str = "./slurm_config.toml"):
-        raw_toml = cls.get_raw_toml(path)
-        parsed_toml = cls.set_correct_config(raw_toml)
-        cls.set_config_values(parsed_toml)
+    def __init__(self, path: Path = default_slurm_path) -> None:
+        raw_toml = self.get_raw_toml(path)
+        self.raw_toml = self.flatten(raw_toml)
 
-    @classmethod
-    def set_correct_config(cls, config_table: dict) -> dict:
-        """Adjust invalid or mutually exclusive configuration values.
-
-        Args:
-            config_table: Loaded .toml data.
+    def validate(self) -> SlurmConfigModel:
+        """Validate configuration using Pydantic class.
 
         Returns:
-            Adjusted .toml data.
+            CosipyConfigModel: Validated configuration.
         """
-        if config_table["OVERRIDES"]["memory_per_process"]:
-            memory = (
-                config_table["OVERRIDES"]["memory_per_process"]
-                * config_table["MEMORY"]["cores"]
-            )
-            config_table["MEMORY"]["memory"] = f"{memory}GB"
-
-        return config_table
+        return SlurmConfigModel(**self.raw_toml)
 
 
-def main():
-    cfg = Config()
-    if cfg.slurm_use:
-        SlurmConfig()
+def main() -> tuple[CosipyConfigModel, Optional[SlurmConfigModel]]:
+    args = get_user_arguments()
+    cfg = Config(args.config_path).validate()
+    slurm_cfg = SlurmConfig(args.slurm_path).validate() if cfg.slurm_use else None
+    return cfg, slurm_cfg
 
 
 if __name__ == "__main__":
     main()
 else:
-    main()
+    main_config, slurm_config = main()
diff --git a/cosipy/constants.py b/cosipy/constants.py
index 223c1eb..1bc9fb8 100644
--- a/cosipy/constants.py
+++ b/cosipy/constants.py
@@ -2,12 +2,98 @@
 from pathlib import Path
 from typing import Literal
 
-from cosipy.config import Config, TomlLoader, get_user_arguments
+from pydantic import BaseModel, model_validator
+from typing_extensions import Self
+
+from cosipy.config import (TomlLoader, default_constants_path,
+                           get_user_arguments, main_config)
 
 if sys.version_info >= (3, 11):
     import tomllib
 else:
-    import tomli as tomllib  # backwards compatibility
+    pass  # backwards compatibility
+
+
+class ConstantsModel(BaseModel):
+    dt: int
+    max_layers: int
+    z: float
+    zlt1: float
+    zlt2: float
+    t_star_K: int
+    t_star_cutoff: float
+    t_star_wet: int
+    t_star_dry: int
+    sfc_temperature_method: str
+    stability_correction: str
+    albedo_method: str
+    densification_method: str
+    penetrating_method: str
+    roughness_method: str
+    saturation_water_vapour_method: str
+    initial_snowheight_constant: float
+    initial_snow_layer_heights: float
+    initial_glacier_height: float
+    initial_glacier_layer_heights: float
+    initial_top_density_snowpack: float
+    initial_bottom_density_snowpack: float
+    temperature_bottom: float
+    const_init_temp: float
+    center_snow_transfer_function: float
+    spread_snow_transfer_function: float
+    mult_factor_RRR: float
+    minimum_snow_layer_height: float
+    minimum_snowfall: float
+    remesh_method: str
+    thermal_conductivity_method: Literal["bulk", "empirical"]
+    first_layer_height: float
+    layer_stretching: float
+    density_threshold_merging: float
+    temperature_threshold_merging: float
+    albedo_fresh_snow: float
+    albedo_firn: float
+    albedo_ice: float
+    albedo_mod_snow_aging: float
+    albedo_mod_snow_depth: float
+    roughness_fresh_snow: float
+    roughness_ice: float
+    roughness_firn: float
+    aging_factor_roughness: float
+    snow_ice_threshold: float
+    lat_heat_melting: float
+    lat_heat_vaporize: float
+    lat_heat_sublimation: float
+    spec_heat_air: float
+    spec_heat_ice: float
+    spec_heat_water: float
+    k_i: float
+    k_w: float
+    k_a: float
+    water_density: float
+    ice_density: float
+    air_density: float
+    sigma: float
+    zero_temperature: float
+    surface_emission_coeff: float
+    merge_max: int
+    constant_density: float
+
+    @model_validator(mode="after")
+    def adjust_config(self) -> Self:
+        """Adjust invalid or mutually exclusive configuration values.
+
+        Args:
+            config_table: Loaded .toml data.
+
+        Returns:
+            Adjusted .toml data.
+        """
+        # WRF_X_CSPY: for efficiency and consistency
+        if main_config.WRF_X_CSPY:
+            self.albedo_method = "Oerlemans98"
+            self.stability_correction = "MO"
+
+        return self
 
 
 class Constants(TomlLoader):
@@ -17,42 +103,29 @@ class Constants(TomlLoader):
     file.
     """
 
-    def __init__(self):
-        self.args = get_user_arguments()
-        self.load(self.args.constants_path)
-
-    @classmethod
-    def load(cls, path: str = "./constants.toml"):
-        raw_toml = cls.get_raw_toml(path)
-        parsed_toml = cls.set_correct_config(raw_toml)
-        cls.set_config_values(parsed_toml)
+    def __init__(self, path: Path = default_constants_path) -> None:
+        raw_toml = self.get_raw_toml(path)
+        self.raw_toml = self.flatten(raw_toml)
 
     @classmethod
-    def set_correct_config(cls, config_table: dict) -> dict:
-        """Adjust invalid or mutually exclusive configuration values.
+    def set_config(cls, config: dict) -> None:
+        for key, value in config.items():
+            setattr(cls, key, value)
 
-        Args:
-            config_table: Loaded .toml data.
+    def validate(self) -> ConstantsModel:
+        """Validate configuration using Pydantic class.
 
         Returns:
-            Adjusted .toml data.
+            CosipyConfigModel: Validated configuration.
         """
-        # WRF_X_CSPY: for efficiency and consistency
-        if Config.WRF_X_CSPY:
-            config_table["PARAMETERIZATIONS"]["albedo_method"] = "Oerlemans98"
-            config_table["PARAMETERIZATIONS"]["stability_correction"] = "MO"
-            # config_table["PARAMETERIZATIONS"][
-            #     "sfc_temperature_method"
-            # ] = "Secant"
-
-        return config_table
-
+        return ConstantsModel(**self.raw_toml)
 
 def main():
-    Constants()
+    args = get_user_arguments()
+    return Constants(args.constants_path).validate()
 
 
 if __name__ == "__main__":
     main()
 else:
-    main()
+    constants_config = main()
diff --git a/cosipy/cpkernel/cosipy_core.py b/cosipy/cpkernel/cosipy_core.py
index b4aa2fb..9a51977 100644
--- a/cosipy/cpkernel/cosipy_core.py
+++ b/cosipy/cpkernel/cosipy_core.py
@@ -1,8 +1,8 @@
 import numpy as np
 import pandas as pd
 
-from cosipy.config import Config
-from cosipy.constants import Constants
+from cosipy.config import main_config
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.init import init_snowpack, load_snowpack
 from cosipy.cpkernel.io import IOClass
 from cosipy.modules.albedo import updateAlbedo
@@ -25,7 +25,7 @@ def init_nan_array_1d(nt: int) -> np.ndarray:
     Returns:
         NaN array.
     """
-    if not Config.WRF_X_CSPY:
+    if not main_config.WRF_X_CSPY:
         x = np.full(nt, np.nan)
     else:
         x = None
@@ -43,7 +43,7 @@ def init_nan_array_2d(nt: int, max_layers: int) -> np.ndarray:
     Returns:
         2D NaN array.
     """
-    if not Config.WRF_X_CSPY and Config.full_field:
+    if not main_config.WRF_X_CSPY and main_config.full_field:
         x = np.full((nt, max_layers), np.nan)
     else:
         x = None
@@ -70,24 +70,24 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
     """
 
     # Declare locally for faster lookup
-    dt = Constants.dt
-    max_layers = Constants.max_layers
-    z = Constants.z
-    mult_factor_RRR = Constants.mult_factor_RRR
-    densification_method = Constants.densification_method
-    ice_density = Constants.ice_density
-    water_density = Constants.water_density
-    minimum_snowfall = Constants.minimum_snowfall
-    zero_temperature = Constants.zero_temperature
-    lat_heat_sublimation = Constants.lat_heat_sublimation
-    lat_heat_melting = Constants.lat_heat_melting
-    lat_heat_vaporize = Constants.lat_heat_vaporize
-    center_snow_transfer_function = Constants.center_snow_transfer_function
-    spread_snow_transfer_function = Constants.spread_snow_transfer_function
-    constant_density = Constants.constant_density
-    albedo_fresh_snow = Constants.albedo_fresh_snow
-    albedo_firn = Constants.albedo_firn
-    WRF_X_CSPY = Config.WRF_X_CSPY
+    dt = cc.dt
+    max_layers = cc.max_layers
+    z = cc.z
+    mult_factor_RRR = cc.mult_factor_RRR
+    densification_method = cc.densification_method
+    ice_density = cc.ice_density
+    water_density = cc.water_density
+    minimum_snowfall = cc.minimum_snowfall
+    zero_temperature = cc.zero_temperature
+    lat_heat_sublimation = cc.lat_heat_sublimation
+    lat_heat_melting = cc.lat_heat_melting
+    lat_heat_vaporize = cc.lat_heat_vaporize
+    center_snow_transfer_function = cc.center_snow_transfer_function
+    spread_snow_transfer_function = cc.spread_snow_transfer_function
+    constant_density = cc.constant_density
+    albedo_fresh_snow = cc.albedo_fresh_snow
+    albedo_firn = cc.albedo_firn
+    WRF_X_CSPY = main_config.WRF_X_CSPY
 
     # Replace values from constants.py if coupled
     # TODO: This only affects the current module scope instead of global.
@@ -188,7 +188,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
         RRR = DATA.RRR.values * mult_factor_RRR
 
     # Use RRR rather than snowfall?
-    if Config.force_use_TP:
+    if main_config.force_use_TP:
         SNOWF = None
 
     LWin = np.array(nt * [None])
@@ -203,7 +203,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
         N = DATA.N.values
 
     # Use N rather than LWin
-    if Config.force_use_N:
+    if main_config.force_use_N:
         LWin = None
 
     SLOPE = 0.
@@ -219,7 +219,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
     if WRF_X_CSPY:
         albedo_snow = albedo_firn
 
-    if Config.stake_evaluation:
+    if main_config.stake_evaluation:
         # Create pandas dataframe for stake evaluation
         _df = pd.DataFrame(index=stake_data.index, columns=['mb','snowheight'], dtype='float')
 
@@ -432,7 +432,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
             _MOL[t] = MOL
             _new_snow_height[t], _new_snow_timestamp[t], _old_snow_timestamp[t] = GRID.get_fresh_snow_props()
 
-            if Config.full_field:
+            if main_config.full_field:
                 if GRID.get_number_layers()>max_layers:
                     raise ValueError('Maximum number of layers reached')
                 _LAYER_HEIGHT[t, 0:GRID.get_number_layers()] = GRID.get_height()
@@ -483,7 +483,7 @@ def cosipy_core(DATA, indY, indX, GRID_RESTART=None, stake_names=None, stake_dat
                     _LAYER_RHO,_LAYER_T,_LAYER_LWC,None,None,_LAYER_ICE_FRACTION,None,None,None,None,None)
 
     if not WRF_X_CSPY:
-        if Config.stake_evaluation:
+        if main_config.stake_evaluation:
             # Evaluate stakes
             _stat = evaluate(stake_names, stake_data, _df)
 
diff --git a/cosipy/cpkernel/grid.py b/cosipy/cpkernel/grid.py
index 94072a1..9665df7 100644
--- a/cosipy/cpkernel/grid.py
+++ b/cosipy/cpkernel/grid.py
@@ -5,7 +5,7 @@
 from numba import float64, intp, optional, typed, types
 from numba.experimental import jitclass
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.node import Node
 
 node_type = Node.class_type.instance_type
@@ -23,18 +23,18 @@
 spec["grid"] = types.ListType(node_type)
 
 # only required for njitted functions
-snow_ice_threshold = Constants.snow_ice_threshold
-first_layer_height = Constants.first_layer_height
-layer_stretching = Constants.layer_stretching
-temperature_threshold_merging = Constants.temperature_threshold_merging
-density_threshold_merging = Constants.density_threshold_merging
-merge_max = Constants.merge_max
-minimum_snowfall = Constants.minimum_snowfall
-minimum_snow_layer_height = Constants.minimum_snow_layer_height
-remesh_method = Constants.remesh_method
-ice_density = Constants.ice_density
-water_density = Constants.water_density
-albedo_method = Constants.albedo_method
+snow_ice_threshold = cc.snow_ice_threshold
+first_layer_height = cc.first_layer_height
+layer_stretching = cc.layer_stretching
+temperature_threshold_merging = cc.temperature_threshold_merging
+density_threshold_merging = cc.density_threshold_merging
+merge_max = cc.merge_max
+minimum_snowfall = cc.minimum_snowfall
+minimum_snow_layer_height = cc.minimum_snow_layer_height
+remesh_method = cc.remesh_method
+ice_density = cc.ice_density
+water_density = cc.water_density
+albedo_method = cc.albedo_method
 
 
 @jitclass(spec)
diff --git a/cosipy/cpkernel/init.py b/cosipy/cpkernel/init.py
index 264b5ec..83e35d2 100644
--- a/cosipy/cpkernel/init.py
+++ b/cosipy/cpkernel/init.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numba import njit
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.grid import Grid
 
 
@@ -16,13 +16,13 @@ def init_snowpack(DATA):
     """
 
     # Declare locally for faster lookup
-    initial_snowheight_constant = Constants.initial_snowheight_constant
-    initial_glacier_height = Constants.initial_glacier_height
-    initial_glacier_layer_heights = Constants.initial_glacier_layer_heights
-    temperature_bottom = Constants.temperature_bottom
-    initial_top_density_snowpack = Constants.initial_top_density_snowpack
-    initial_bottom_density_snowpack = Constants.initial_bottom_density_snowpack
-    ice_density = Constants.ice_density
+    initial_snowheight_constant = cc.initial_snowheight_constant
+    initial_glacier_height = cc.initial_glacier_height
+    initial_glacier_layer_heights = cc.initial_glacier_layer_heights
+    temperature_bottom = cc.temperature_bottom
+    initial_top_density_snowpack = cc.initial_top_density_snowpack
+    initial_bottom_density_snowpack = cc.initial_bottom_density_snowpack
+    ice_density = cc.ice_density
 
     # Check for WRF data
     if "SNOWHEIGHT" in DATA:
@@ -31,7 +31,7 @@ def init_snowpack(DATA):
             initial_snowheight = 0.0
     else:
         initial_snowheight = initial_snowheight_constant
-    temperature_top = np.minimum(DATA.T2.values[0], Constants.zero_temperature)
+    temperature_top = np.minimum(DATA.T2.values[0], cc.zero_temperature)
 
     # Do the vertical interpolation
     layer_heights = []
diff --git a/cosipy/cpkernel/io.py b/cosipy/cpkernel/io.py
index ea9343f..1cf8a41 100644
--- a/cosipy/cpkernel/io.py
+++ b/cosipy/cpkernel/io.py
@@ -10,8 +10,8 @@
 import numpy as np
 import xarray as xr
 
-from cosipy.config import Config
-from cosipy.constants import Constants
+from cosipy.config import main_config
+from cosipy.constants import constants_config as cc
 
 
 class IOClass:
@@ -23,18 +23,17 @@ def __init__(self, data: xr.Dataset | None = None):
             RESTART (xarray.Dataset or None): Restart file data.
             RESULT (xarray.Dataset or None): Model result data.
         """
-
-        self.atm = self.get_output_variables(Config.output_atm)
-        self.internal = self.get_output_variables(Config.output_internal)
-        self.full = self.get_output_variables(Config.output_full)
+        self.atm = self.get_output_variables(main_config.output_atm)
+        self.internal = self.get_output_variables(main_config.output_internal)
+        self.full = self.get_output_variables(main_config.output_full)
 
         # Initialize data
-        self.DATA = DATA
+        self.DATA = data
         self.RESTART = None
         self.RESULT = None
 
         # If local IO class is initialized we need to get the dimensions of the dataset
-        if DATA is not None:
+        if self.DATA is not None:
             self.time = self.DATA.sizes["time"]
 
     def get_output_variables(self, variables: str) -> list:
@@ -42,25 +41,24 @@ def get_output_variables(self, variables: str) -> list:
         # Sets are unordered -> same config, different checksums
         if not variables:
             return []
-        else:
-            return [item for item in variables.split(",")]
+        return list(variables.split(","))
 
-    def init_atm_attribute(self, name: str):
+    def init_atm_attribute(self, name: str) -> None:
         """Initialise empty atm attribute."""
         if name in self.atm:
             setattr(self, name, self.create_nan_array())
 
-    def init_internal_attribute(self, name: str):
+    def init_internal_attribute(self, name: str) -> None:
         """Initialise empty internal attribute."""
         if name in self.internal:
             setattr(self, name, self.create_nan_array())
 
-    def init_full_field_attribute(self, name: str, max_layers: int):
+    def init_full_field_attribute(self, name: str, max_layers: int) -> None:
         """Initialise empty layer attribute."""
         if name in self.full:
             setattr(self, f"LAYER_{name}", self.create_3d_nan_array(max_layers))
 
-    def set_atm_attribute(self, name: str, value: np.ndarray, x: int, y: int):
+    def set_atm_attribute(self, name: str, value: np.ndarray, x: int, y: int) -> None:
         """Set atm attribute if it is a desired output variable.
 
         Args:
@@ -72,7 +70,9 @@ def set_atm_attribute(self, name: str, value: np.ndarray, x: int, y: int):
         if name in self.atm:
             getattr(self, name)[:, y, x] = value
 
-    def set_internal_attribute(self, name: str, value: np.ndarray, x: int, y: int):
+    def set_internal_attribute(
+        self, name: str, value: np.ndarray, x: int, y: int
+    ) -> None:
         """Set internal attribute if it is a desired output variable.
 
         Args:
@@ -84,7 +84,9 @@ def set_internal_attribute(self, name: str, value: np.ndarray, x: int, y: int):
         if name in self.internal:
             getattr(self, name)[:, y, x] = value
 
-    def set_full_field_attribute(self, name: str, value: np.ndarray, x: int, y: int):
+    def set_full_field_attribute(
+        self, name: str, value: np.ndarray, x: int, y: int
+    ) -> None:
         """Set layer attribute if it is a desired output variable.
 
         Args:
@@ -97,8 +99,12 @@ def set_full_field_attribute(self, name: str, value: np.ndarray, x: int, y: int)
             getattr(self, f"LAYER_{name}")[:, y, x, :] = value
 
     def get_datetime(
-        self, timestamp: str, use_np: bool = True, fmt: str = "%Y-%m-%dT%H:%M"
-    ):
+        self,
+        timestamp: str | datetime | np.datetime64,
+        *,
+        use_np: bool = True,
+        fmt: str = "%Y-%m-%dT%H:%M",
+    ) -> datetime | np.datetime64:
         """Get datetime object from a string.
 
         Args:
@@ -113,10 +119,8 @@ def get_datetime(
         if isinstance(timestamp, str):
             if use_np:
                 return np.datetime64(timestamp)
-            else:
-                return datetime.strptime(timestamp, fmt)
-        else:
-            return timestamp
+            return datetime.strptime(timestamp, fmt)
+        return timestamp
 
     def create_data_file(self) -> xr.Dataset:
         """Create the input data and read the restart file if necessary.
@@ -124,15 +128,16 @@ def create_data_file(self) -> xr.Dataset:
         Returns:
             Model input data.
         """
-
-        if Config.restart:
+        if not self.DATA:
+            self.DATA = self.create_data_file()
+        if main_config.restart:
             print(f"{'-'*62}\n\tRESTART FROM PREVIOUS STATE\n{'-'*62}\n")
 
             # Load the restart file
-            time_start = Config.time_start
-            time_end = Config.time_end
-            start_timestamp = self.get_datetime(time_start)
-            end_timestamp = self.get_datetime(time_end)
+            time_start = main_config.time_start
+            time_end = main_config.time_end
+            start_timestamp = datetime.strptime(time_start, "%Y-%m-%dT%H:%M%z")
+            end_timestamp = datetime.strptime(time_end, "%Y-%m-%dT%H:%M%z")
             timestamp = start_timestamp.strftime("%Y-%m-%dT%H-%M")
             restart_path = Path(Config.data_path) / "restart" / f"restart_{timestamp}.nc"
             if not restart_path.is_file():
@@ -142,44 +147,41 @@ def create_data_file(self) -> xr.Dataset:
                 msg = f"Start date {time_start} equals end date {time_end}"
                 raise IndexError(msg)
             try:
-                if not restart_path.is_file():
-                    raise FileNotFoundError
-                elif start_timestamp == end_timestamp:
-                    raise IndexError
-                else:
-                    self.GRID_RESTART = xr.open_dataset(restart_path)
-                    """Get time of the last calculation and add one time
-                    step. GRID_RESTART.time is an array of np.datetime64
-                    objects.
-                    """
-                    self.restart_date = self.GRID_RESTART.time.values + np.timedelta64(
-                        Constants.dt, "s"
-                    )
-                    # Read data from the last date to the end of the data file
-                    self.init_data_dataset()
+                self.GRID_RESTART = xr.open_dataset(restart_path)
+                """Get time of the last calculation and add one time step.
+
+                GRID_RESTART.time is an array of np.datetime64 objects.
+                """
+                self.restart_date = self.GRID_RESTART.time.to_numpy() + np.timedelta64(
+                    cc.dt, "s"
+                )
+                # Read data from the last date to the end of the data file
+                self.init_data_dataset()
             except FileNotFoundError:
-                raise SystemExit(f"No restart file available for the given date: {timestamp}")
+                msg = f"No restart file available for the given date: {timestamp}"
+                raise SystemExit(msg) from None
             except IndexError:
-                raise SystemExit(f"Start date {time_start} equals end date {time_end}\n")
+                msg = f"Start date {time_start} equals end date {time_end}\n"
+                raise SystemExit(msg) from None
         else:
             # If no restart, read data according to the dates defined in config file
             self.restart_date = None
             self.init_data_dataset()
 
-        if Config.tile:  # Tile the data if desired
-            if Config.WRF:
+        if main_config.tile:  # Tile the data if desired
+            if main_config.WRF:
                 self.DATA = self.DATA.isel(
-                    south_north=slice(Config.ystart, Config.yend),
-                    west_east=slice(Config.xstart, Config.xend),
+                    south_north=slice(main_config.ystart, main_config.yend),
+                    west_east=slice(main_config.xstart, main_config.xend),
                 )
             else:
                 self.DATA = self.DATA.isel(
-                    lat=slice(Config.ystart, Config.yend),
-                    lon=slice(Config.xstart, Config.xend),
+                    lat=slice(main_config.ystart, main_config.yend),
+                    lon=slice(main_config.xstart, main_config.xend),
                 )
 
-        self.ny = self.DATA.sizes[Config.northing]
-        self.nx = self.DATA.sizes[Config.easting]
+        self.ny = self.DATA.sizes[main_config.northing]
+        self.nx = self.DATA.sizes[main_config.easting]
         self.time = self.DATA.sizes["time"]
 
         return self.DATA
@@ -194,7 +196,7 @@ def create_restart_file(self) -> xr.Dataset:
         self.init_restart_dataset()
         return self.RESTART
 
-    def create_grid_restart(self):
+    def create_grid_restart(self) -> xr.Dataset:
         """Create and initialise the GRID_RESTART structure.
 
         This contains the layer state of the last time step, which is
@@ -208,7 +210,6 @@ def create_nan_array(self) -> np.ndarray:
         Returns:
             Filled array with time and 2D spatial coordinates.
         """
-
         return np.full((self.time, self.ny, self.nx), np.nan)
 
     def create_3d_nan_array(self, max_layers: int) -> np.ndarray:
@@ -220,7 +221,6 @@ def create_3d_nan_array(self, max_layers: int) -> np.ndarray:
         Returns:
             Filled array with time and 3D spatial coordinates.
         """
-
         return np.full((self.time, self.ny, self.nx, max_layers), np.nan)
 
     def check_field(self, field, _max, _min) -> bool:
@@ -232,13 +232,12 @@ def check_field(self, field, _max, _min) -> bool:
                 f"MIN: {np.nanmin(field):.2f}\n"
             )
             return False
-        else:
-            return True
+        return True
 
     def check_input_data(self) -> bool:
         """Check the input data is within valid bounds."""
         print(f"{'-'*62}\nChecking input data ....\n")
-        
+
         data_bounds = {
             "T2": (313.16, 243.16),
             "RH2": (100.0, 0.0),
@@ -246,24 +245,25 @@ def check_input_data(self) -> bool:
             "U2": (50.0, 0.0),
             "RRR": (20.0, 0.0),
             "N": (1.0, 0.0),
-            "PRES": (1080.0, 400.0),
+            "PRESS": (1080.0, 400.0),
             "LWin": (400.0, 200.0),
             "SNOWFALL": (0.1, 0.0),
             "SLOPE": (0.0, 90.0),
         }
 
         for key, bounds in data_bounds.items():
-            if key in self.DATA:
-                if self.check_field(self.DATA[key], bounds[0], bounds[1]):
-                    print(f"{key} ... ok")
+            if key in self.DATA and self.check_field(
+                self.DATA[key], bounds[0], bounds[1]
+            ):
+                print(f"{key} ... ok")
 
         return True
 
-    def init_data_dataset(self):
+    def init_data_dataset(self) -> None:
         """Read and store the input netCDF data.
 
         The input data should contain the following variables:
-            :PRES: Air pressure [hPa].
+            :PRESS: Air pressure [hPa].
             :N: Cloud cover fraction [-].
             :RH2: 2m relative humidity [%].
             :RRR: Precipitation per time step [mm].
@@ -273,16 +273,16 @@ def init_data_dataset(self):
             :U2: Wind speed (magnitude) [|m s^-1|].
             :HGT: Elevation [m].
         """
-        input_path = Path(Config.data_path) / "input" / Config.input_netcdf
+        input_path = Path(main_config.data_path) / "input" / main_config.input_netcdf
         try:
             self.DATA = xr.open_dataset(input_path)
-        except FileNotFoundError:
-            raise SystemExit(f"Input file not found at: {input_path}")
-
+        except FileNotFoundError as e:
+            msg = f"Input file not found at: {input_path}"
+            raise SystemExit(msg) from e
 
         self.DATA["time"] = np.sort(self.DATA["time"].values)
-        minimum_time = str(self.DATA.time.values[0])[0:16]
-        maximum_time = str(self.DATA.time.values[-1])[0:16]
+        minimum_time = str(self.DATA.time.to_numpy()[0])[0:16]
+        maximum_time = str(self.DATA.time.to_numpy()[-1])[0:16]
         start_interval = self.get_datetime(minimum_time)
         end_interval = self.get_datetime(maximum_time)
         time_steps = str(self.DATA.sizes["time"])
@@ -291,38 +291,37 @@ def init_data_dataset(self):
             f"until {maximum_time}. Time steps: {time_steps}\n\n"
         )
 
-        time_start = Config.time_start  # avoid repeat calls
-        time_end = Config.time_end
+        time_start = main_config.time_start  # avoid repeat calls
+        time_end = main_config.time_end
         start_time = self.get_datetime(time_start)
         end_time = self.get_datetime(time_end)
 
         if (start_time > end_interval) or (end_time < start_interval):
-            raise IndexError("Selected period not available in input data.\n")
+            msg = "Selected period not available in input data.\n"
+            raise IndexError(msg)
         if start_time < start_interval:
             warnings.warn(
-                "\nWARNING! Selected startpoint before first timestep of input data\n",    
+                "\nWARNING! Selected startpoint before first timestep of input data\n",
             )
         if end_time > end_interval:
             warnings.warn(
-                "\nWARNING! Selected endpoint after last timestep of input data\n",    
+                "\nWARNING! Selected endpoint after last timestep of input data\n",
             )
 
         if self.restart_date is None:  # Check if restart option is set
-            print(
-                f"{'-'*62}\n\tIntegration from {time_start} to {time_end}\n{'-'*62}\n"
-            )
+            print(f"{'-'*62}\n\tIntegration from {time_start} to {time_end}\n{'-'*62}\n")
             self.DATA = self.DATA.sel(time=slice(time_start, time_end))
+        # There is nothing to do if the dates are equal
+        elif self.restart_date == end_time:
+            msg = "Start date equals end date ... no new data ... EXIT"
+            raise SystemExit(msg)
         else:
-            # There is nothing to do if the dates are equal
-            if self.restart_date == end_time:
-                raise SystemExit("Start date equals end date ... no new data ... EXIT")
-            else:
-                # otherwise, run the model from the restart date to the defined end date
-                print(
-                    f"Starting from {self.restart_date} (from restart file) "
-                    f"to {time_end} (from config file)\n"
-                )
-                self.DATA = self.DATA.sel(time=slice(self.restart_date, time_end))
+            # otherwise, run the model from the restart date to the defined end date
+            print(
+                f"Starting from {self.restart_date} (from restart file) "
+                f"to {time_end} (from config file)\n"
+            )
+            self.DATA = self.DATA.sel(time=slice(self.restart_date, time_end))
 
         self.check_input_data()
         print(f"\nGlacier gridpoints: {np.nansum(self.DATA.MASK >= 1)} \n\n")
@@ -344,7 +343,7 @@ def get_input_metadata(self) -> tuple:
             "T2": ("K", "Air temperature at 2 m"),
             "RH2": ("%", "Relative humidity at 2 m"),
             "U2": ("m s\u207b\xb9", "Wind velocity at 2 m"),
-            "PRES": ("hPa", "Atmospheric pressure"),
+            "PRESS": ("hPa", "Atmospheric pressure"),
             "G": ("W m\u207b\xb2", "Incoming shortwave radiation"),
             "RRR": ("mm", "Total precipitation"),
             "SNOWFALL": ("m", "Snowfall"),
@@ -354,14 +353,14 @@ def get_input_metadata(self) -> tuple:
 
         return metadata_spatial, metadata_spatiotemporal
 
-    def get_full_field_metadata(self) -> dict:
+    def get_full_field_metadata(self) -> dict[str, tuple[str, str]]:
         """Get layer variable names and units.
 
         Returns:
             Metadata for full field layer variables, including netCDF
             keyname, unit, and long name.
         """
-        metadata = {
+        return {
             "LAYER_HEIGHT": ("m", "Layer height"),
             "LAYER_RHO": ("kg m^-3", "Layer density"),
             "LAYER_T": ("K", "Layer temperature"),
@@ -374,10 +373,8 @@ def get_full_field_metadata(self) -> dict:
             "LAYER_REFREEZE": ("m w.e.", "Refreezing"),
         }
 
-        return metadata
-
-    def get_restart_metadata(self) -> dict:
 
+    def get_restart_metadata(self) -> dict[str, tuple[str, str]]:
         field_metadata = self.get_full_field_metadata()
         restart_metadata = {
             "new_snow_height": ("m .w.e", "New snow height"),
@@ -388,11 +385,10 @@ def get_restart_metadata(self) -> dict:
             "NEWSNOWTIMESTAMP": ("s", "New snow timestamp"),
             "OLDSNOWTIMESTAMP": ("s", "Old snow timestamp"),
         }
-        metadata = restart_metadata | field_metadata
+        return restart_metadata | field_metadata
 
-        return metadata
 
-    def get_result_metadata(self) -> dict:
+    def get_result_metadata(self) -> dict[str, tuple[str, str]]:
         """Get all variable names and units.
 
         Returns:
@@ -431,11 +427,10 @@ def get_result_metadata(self) -> dict:
             "MOL": ("m", "Monin Obukhov length"),
         }
         # metadata_result overwrites items in previous union!
-        metadata = (
+        return (
             metadata_spatial | metadata_spatiotemporal | metadata_full | metadata_result
         )
 
-        return metadata
 
     def init_result_dataset(self) -> xr.Dataset:
         """Create the final dataset to aggregate and store the results.
@@ -447,7 +442,6 @@ def init_result_dataset(self) -> xr.Dataset:
         Returns:
             One-dimensional structure with the model output.
         """
-
         # Coordinates
         self.RESULT = xr.Dataset()
         self.RESULT.coords["time"] = self.DATA.coords["time"]
@@ -455,78 +449,100 @@ def init_result_dataset(self) -> xr.Dataset:
         self.RESULT.coords["lon"] = self.DATA.coords["lon"]
 
         # Global attributes from config.py
-        self.RESULT.attrs["Start_from_restart_file"] = str(Config.restart)
-        self.RESULT.attrs["Stake_evaluation"] = str(Config.stake_evaluation)
-        self.RESULT.attrs["WRF_simulation"] = str(Config.WRF)
-        self.RESULT.attrs["Compression_level"] = Config.compression_level
-        self.RESULT.attrs["Slurm_use"] = str(Config.slurm_use)
-        self.RESULT.attrs["Full_field"] = str(Config.full_field)
-        self.RESULT.attrs["Force_use_TP"] = str(Config.force_use_TP)
-        self.RESULT.attrs["Force_use_N"] = str(Config.force_use_N)
-        self.RESULT.attrs["Tile_of_glacier_of_interest"] = str(Config.tile)
+        self.RESULT.attrs["Start_from_restart_file"] = str(main_config.restart)
+        self.RESULT.attrs["Stake_evaluation"] = str(main_config.stake_evaluation)
+        self.RESULT.attrs["WRF_simulation"] = str(main_config.WRF)
+        self.RESULT.attrs["Compression_level"] = main_config.compression_level
+        self.RESULT.attrs["Slurm_use"] = str(main_config.slurm_use)
+        self.RESULT.attrs["Full_field"] = str(main_config.full_field)
+        self.RESULT.attrs["Force_use_TP"] = str(main_config.force_use_TP)
+        self.RESULT.attrs["Force_use_N"] = str(main_config.force_use_N)
+        self.RESULT.attrs["Tile_of_glacier_of_interest"] = str(main_config.tile)
 
         # Global attributes from constants.py
-        self.RESULT.attrs["Time_step_input_file_seconds"] = Constants.dt
-        self.RESULT.attrs["Max_layers"] = Constants.max_layers
-        self.RESULT.attrs["Z_measurement_height"] = Constants.z
-        self.RESULT.attrs["Stability_correction"] = Constants.stability_correction
-        self.RESULT.attrs["Albedo_method"] = Constants.albedo_method
-        self.RESULT.attrs["Densification_method"] = Constants.densification_method
-        self.RESULT.attrs["Penetrating_method"] = Constants.penetrating_method
-        self.RESULT.attrs["Roughness_method"] = Constants.roughness_method
-        self.RESULT.attrs["Saturation_water_vapour_method"] = Constants.saturation_water_vapour_method
-
-        self.RESULT.attrs["Initial_snowheight"] = Constants.initial_snowheight_constant
-        self.RESULT.attrs["Initial_snow_layer_heights"] = Constants.initial_snow_layer_heights
-        self.RESULT.attrs["Initial_glacier_height"] = Constants.initial_glacier_height
-        self.RESULT.attrs["Initial_glacier_layer_heights"] = Constants.initial_glacier_layer_heights
-        self.RESULT.attrs["Initial_top_density_snowpack"] = Constants.initial_top_density_snowpack
-        self.RESULT.attrs["Initial_bottom_density_snowpack"] = Constants.initial_bottom_density_snowpack
-        self.RESULT.attrs["Temperature_bottom"] = Constants.temperature_bottom
-        self.RESULT.attrs["Const_init_temp"] = Constants.const_init_temp
-
-        self.RESULT.attrs["Center_snow_transfer_function"] = Constants.center_snow_transfer_function
-        self.RESULT.attrs["Spread_snow_transfer_function"] = Constants.spread_snow_transfer_function
-        self.RESULT.attrs["Multiplication_factor_for_RRR_or_SNOWFALL"] = Constants.mult_factor_RRR
-        self.RESULT.attrs["Minimum_snow_layer_height"] = Constants.minimum_snow_layer_height
-        self.RESULT.attrs["Minimum_snowfall"] = Constants.minimum_snowfall
-
-        self.RESULT.attrs["Remesh_method"] = Constants.remesh_method
-        self.RESULT.attrs["First_layer_height_log_profile"] = Constants.first_layer_height
-        self.RESULT.attrs["Layer_stretching_log_profile"] = Constants.layer_stretching
-
-        self.RESULT.attrs["Merge_max"] = Constants.merge_max
-        self.RESULT.attrs["Layer_stretching_log_profile"] = Constants.layer_stretching
-        self.RESULT.attrs["Density_threshold_merging"] = Constants.density_threshold_merging
-        self.RESULT.attrs["Temperature_threshold_merging"] = Constants.temperature_threshold_merging
-
-        self.RESULT.attrs["Density_fresh_snow"] = Constants.constant_density
-        self.RESULT.attrs["Albedo_fresh_snow"] = Constants.albedo_fresh_snow
-        self.RESULT.attrs["Albedo_firn"] = Constants.albedo_firn
-        self.RESULT.attrs["Albedo_ice"] = Constants.albedo_ice
-        self.RESULT.attrs["Albedo_mod_snow_aging"] = Constants.albedo_mod_snow_aging
-        self.RESULT.attrs["Albedo_mod_snow_depth"] = Constants.albedo_mod_snow_depth
-        self.RESULT.attrs["Roughness_fresh_snow"] = Constants.roughness_fresh_snow
-        self.RESULT.attrs["Roughness_ice"] = Constants.roughness_ice
-        self.RESULT.attrs["Roughness_firn"] = Constants.roughness_firn
-        self.RESULT.attrs["Aging_factor_roughness"] = Constants.aging_factor_roughness
-        self.RESULT.attrs["Snow_ice_threshold"] = Constants.snow_ice_threshold
-
-        self.RESULT.attrs["lat_heat_melting"] = Constants.lat_heat_melting
-        self.RESULT.attrs["lat_heat_vaporize"] = Constants.lat_heat_vaporize
-        self.RESULT.attrs["lat_heat_sublimation"] = Constants.lat_heat_sublimation
-        self.RESULT.attrs["spec_heat_air"] = Constants.spec_heat_air
-        self.RESULT.attrs["spec_heat_ice"] = Constants.spec_heat_ice
-        self.RESULT.attrs["spec_heat_water"] = Constants.spec_heat_water
-        self.RESULT.attrs["k_i"] = Constants.k_i
-        self.RESULT.attrs["k_w"] = Constants.k_w
-        self.RESULT.attrs["k_a"] = Constants.k_a
-        self.RESULT.attrs["water_density"] = Constants.water_density
-        self.RESULT.attrs["ice_density"] = Constants.ice_density
-        self.RESULT.attrs["air_density"] = Constants.air_density
-        self.RESULT.attrs["sigma"] = Constants.sigma
-        self.RESULT.attrs["zero_temperature"] = Constants.zero_temperature
-        self.RESULT.attrs["Surface_emission_coeff"] = Constants.surface_emission_coeff
+        self.RESULT.attrs["Time_step_input_file_seconds"] = cc.dt
+        self.RESULT.attrs["Max_layers"] = cc.max_layers
+        self.RESULT.attrs["Z_measurement_height"] = cc.z
+        self.RESULT.attrs["Stability_correction"] = cc.stability_correction
+        self.RESULT.attrs["Albedo_method"] = cc.albedo_method
+        self.RESULT.attrs["Densification_method"] = cc.densification_method
+        self.RESULT.attrs["Penetrating_method"] = cc.penetrating_method
+        self.RESULT.attrs["Roughness_method"] = cc.roughness_method
+        self.RESULT.attrs["Saturation_water_vapour_method"] = (
+            cc.saturation_water_vapour_method
+        )
+
+        self.RESULT.attrs["Initial_snowheight"] = cc.initial_snowheight_constant
+        self.RESULT.attrs["Initial_snow_layer_heights"] = (
+            cc.initial_snow_layer_heights
+        )
+        self.RESULT.attrs["Initial_glacier_height"] = cc.initial_glacier_height
+        self.RESULT.attrs["Initial_glacier_layer_heights"] = (
+            cc.initial_glacier_layer_heights
+        )
+        self.RESULT.attrs["Initial_top_density_snowpack"] = (
+            cc.initial_top_density_snowpack
+        )
+        self.RESULT.attrs["Initial_bottom_density_snowpack"] = (
+            cc.initial_bottom_density_snowpack
+        )
+        self.RESULT.attrs["Temperature_bottom"] = cc.temperature_bottom
+        self.RESULT.attrs["Const_init_temp"] = cc.const_init_temp
+
+        self.RESULT.attrs["Center_snow_transfer_function"] = (
+            cc.center_snow_transfer_function
+        )
+        self.RESULT.attrs["Spread_snow_transfer_function"] = (
+            cc.spread_snow_transfer_function
+        )
+        self.RESULT.attrs["Multiplication_factor_for_RRR_or_SNOWFALL"] = (
+            cc.mult_factor_RRR
+        )
+        self.RESULT.attrs["Minimum_snow_layer_height"] = (
+            cc.minimum_snow_layer_height
+        )
+        self.RESULT.attrs["Minimum_snowfall"] = cc.minimum_snowfall
+
+        self.RESULT.attrs["Remesh_method"] = cc.remesh_method
+        self.RESULT.attrs["First_layer_height_log_profile"] = cc.first_layer_height
+        self.RESULT.attrs["Layer_stretching_log_profile"] = cc.layer_stretching
+
+        self.RESULT.attrs["Merge_max"] = cc.merge_max
+        self.RESULT.attrs["Layer_stretching_log_profile"] = cc.layer_stretching
+        self.RESULT.attrs["Density_threshold_merging"] = (
+            cc.density_threshold_merging
+        )
+        self.RESULT.attrs["Temperature_threshold_merging"] = (
+            cc.temperature_threshold_merging
+        )
+
+        self.RESULT.attrs["Density_fresh_snow"] = cc.constant_density
+        self.RESULT.attrs["Albedo_fresh_snow"] = cc.albedo_fresh_snow
+        self.RESULT.attrs["Albedo_firn"] = cc.albedo_firn
+        self.RESULT.attrs["Albedo_ice"] = cc.albedo_ice
+        self.RESULT.attrs["Albedo_mod_snow_aging"] = cc.albedo_mod_snow_aging
+        self.RESULT.attrs["Albedo_mod_snow_depth"] = cc.albedo_mod_snow_depth
+        self.RESULT.attrs["Roughness_fresh_snow"] = cc.roughness_fresh_snow
+        self.RESULT.attrs["Roughness_ice"] = cc.roughness_ice
+        self.RESULT.attrs["Roughness_firn"] = cc.roughness_firn
+        self.RESULT.attrs["Aging_factor_roughness"] = cc.aging_factor_roughness
+        self.RESULT.attrs["Snow_ice_threshold"] = cc.snow_ice_threshold
+
+        self.RESULT.attrs["lat_heat_melting"] = cc.lat_heat_melting
+        self.RESULT.attrs["lat_heat_vaporize"] = cc.lat_heat_vaporize
+        self.RESULT.attrs["lat_heat_sublimation"] = cc.lat_heat_sublimation
+        self.RESULT.attrs["spec_heat_air"] = cc.spec_heat_air
+        self.RESULT.attrs["spec_heat_ice"] = cc.spec_heat_ice
+        self.RESULT.attrs["spec_heat_water"] = cc.spec_heat_water
+        self.RESULT.attrs["k_i"] = cc.k_i
+        self.RESULT.attrs["k_w"] = cc.k_w
+        self.RESULT.attrs["k_a"] = cc.k_a
+        self.RESULT.attrs["water_density"] = cc.water_density
+        self.RESULT.attrs["ice_density"] = cc.ice_density
+        self.RESULT.attrs["air_density"] = cc.air_density
+        self.RESULT.attrs["sigma"] = cc.sigma
+        self.RESULT.attrs["zero_temperature"] = cc.zero_temperature
+        self.RESULT.attrs["Surface_emission_coeff"] = cc.surface_emission_coeff
 
         # Variables given by the input dataset
         spatial, spatiotemporal = self.get_input_metadata()
@@ -559,18 +575,17 @@ def init_result_dataset(self) -> xr.Dataset:
                 spatiotemporal["N"][0],
             )
 
-        print(f"\nOutput dataset ... ok")
+        print("\nOutput dataset ... ok")
 
         return self.RESULT
 
-    def create_global_result_arrays(self):
+    def create_global_result_arrays(self) -> None:
         """Create the global numpy arrays to store each output variable.
 
         Each global array is filled with local results from the workers.
         The arrays are then assigned to the RESULT dataset and stored to
         disk (see COSIPY.py).
         """
-
         if self.atm:
             for atm_var in self.atm:
                 self.init_atm_attribute(atm_var)
@@ -579,8 +594,8 @@ def create_global_result_arrays(self):
             for internal_var in self.internal:
                 self.init_internal_attribute(internal_var)
 
-        if Config.full_field and self.full:
-            max_layers = Constants.max_layers  # faster lookup
+        if main_config.full_field and self.full:
+            max_layers = cc.max_layers  # faster lookup
             for full_field_var in self.full:
                 self.init_full_field_attribute(full_field_var, max_layers)
 
@@ -624,9 +639,8 @@ def copy_local_to_global(
         local_LAYER_ICE_FRACTION: np.ndarray,
         local_LAYER_IRREDUCIBLE_WATER: np.ndarray,
         local_LAYER_REFREEZE: np.ndarray,
-    ):
+    ) -> None:
         """Copy the local results from workers to global numpy arrays."""
-
         self.set_atm_attribute("RAIN", local_RAIN, x, y)
         self.set_atm_attribute("SNOWFALL", local_SNOWFALL, x, y)
         self.set_atm_attribute("LWin", local_LWin, x, y)
@@ -656,24 +670,21 @@ def copy_local_to_global(
         self.set_internal_attribute("surfM", local_surfM, x, y)
         self.set_internal_attribute("MOL", local_MOL, x, y)
 
-        if Config.full_field:
+        if main_config.full_field:
             self.set_full_field_attribute("HEIGHT", local_LAYER_HEIGHT, x, y)
             self.set_full_field_attribute("RHO", local_LAYER_RHO, x, y)
             self.set_full_field_attribute("T", local_LAYER_T, x, y)
             self.set_full_field_attribute("LWC", local_LAYER_LWC, x, y)
             self.set_full_field_attribute("CC", local_LAYER_CC, x, y)
             self.set_full_field_attribute("POROSITY", local_LAYER_POROSITY, x, y)
-            self.set_full_field_attribute(
-                "ICE_FRACTION", local_LAYER_ICE_FRACTION, x, y
-            )
+            self.set_full_field_attribute("ICE_FRACTION", local_LAYER_ICE_FRACTION, x, y)
             self.set_full_field_attribute(
                 "IRREDUCIBLE_WATER", local_LAYER_IRREDUCIBLE_WATER, x, y
             )
             self.set_full_field_attribute("REFREEZE", local_LAYER_REFREEZE, x, y)
 
-    def write_results_to_file(self):
+    def write_results_to_file(self) -> None:
         """Add the global numpy arrays to the RESULT dataset."""
-
         metadata = self.get_result_metadata()
         if self.atm:
             for atm_var in self.atm:
@@ -695,16 +706,16 @@ def write_results_to_file(self):
                     metadata[internal_var][1],
                 )
 
-        if Config.full_field and self.full:
-                for full_field_var in self.full:
-                    layer_name = f"LAYER_{full_field_var}"
-                    self.add_variable_along_latlonlayertime(
-                        self.RESULT,
-                        getattr(self, layer_name),
-                        layer_name,
-                        metadata[layer_name][0],
-                        metadata[layer_name][1],
-                    )
+        if main_config.full_field and self.full:
+            for full_field_var in self.full:
+                layer_name = f"LAYER_{full_field_var}"
+                self.add_variable_along_latlonlayertime(
+                    self.RESULT,
+                    getattr(self, layer_name),
+                    layer_name,
+                    metadata[layer_name][0],
+                    metadata[layer_name][1],
+                )
 
     def create_empty_restart(self) -> xr.Dataset:
         """Create an empty dataset for the RESTART attribute.
@@ -716,7 +727,7 @@ def create_empty_restart(self) -> xr.Dataset:
         dataset.coords["time"] = self.DATA.coords["time"][-1]
         dataset.coords["lat"] = self.DATA.coords["lat"]
         dataset.coords["lon"] = self.DATA.coords["lon"]
-        dataset.coords["layer"] = np.arange(Constants.max_layers)
+        dataset.coords["layer"] = np.arange(cc.max_layers)
 
         return dataset
 
@@ -731,15 +742,14 @@ def init_restart_dataset(self) -> xr.Dataset:
 
         return self.RESTART
 
-    def create_global_restart_arrays(self):
+    def create_global_restart_arrays(self) -> None:
         """Initialise the global numpy arrays to store layer profiles.
 
         Each global array will be filled with local results from the
         workers. The arrays will then be assigned to the RESTART dataset
         and stored to disk (see COSIPY.py).
         """
-
-        max_layers = Constants.max_layers  # faster lookup
+        max_layers = cc.max_layers  # faster lookup
 
         for name in [
             "NLAYERS",
@@ -762,7 +772,6 @@ def create_local_restart_dataset(self) -> xr.Dataset:
         Returns:
             RESTART dataset initialised with layer profiles.
         """
-
         self.RESTART = self.create_empty_restart()
 
         metadata = self.get_restart_metadata()
@@ -792,7 +801,7 @@ def create_local_restart_dataset(self) -> xr.Dataset:
 
         return self.RESTART
 
-    def copy_local_restart_to_global(self, y: int, x: int, local_restart: xr.Dataset):
+    def copy_local_restart_to_global(self, y: int, x: int, local_restart: xr.Dataset) -> None:
         """Copy local restart data from workers to global numpy arrays.
 
         Args:
@@ -800,7 +809,6 @@ def copy_local_restart_to_global(self, y: int, x: int, local_restart: xr.Dataset
             x: Longitude index.
             local_restart: Local RESTART dataset.
         """
-
         for name in [
             "NLAYERS",
             "NEWSNOWHEIGHT",
@@ -814,9 +822,8 @@ def copy_local_restart_to_global(self, y: int, x: int, local_restart: xr.Dataset
                 local_restart, f"LAYER_{name}"
             )
 
-    def write_restart_to_file(self):
+    def write_restart_to_file(self) -> None:
         """Add global numpy arrays to the RESTART dataset."""
-
         metadata = self.get_restart_metadata()
         for name in [
             "NLAYERS",
@@ -846,27 +853,35 @@ def write_restart_to_file(self):
     @property
     def RAIN(self):
         return self.__RAIN
+
     @property
     def SNOWFALL(self):
         return self.__SNOWFALL
+
     @property
     def LWin(self):
         return self.__LWin
+
     @property
     def LWout(self):
         return self.__LWout
+
     @property
     def H(self):
         return self.__H
+
     @property
     def LE(self):
         return self.__LE
+
     @property
     def B(self):
         return self.__B
+
     @property
     def QRR(self):
         return self.__QRR
+
     @property
     def MB(self):
         return self.__MB
@@ -874,27 +889,35 @@ def MB(self):
     @RAIN.setter
     def RAIN(self, x):
         self.__RAIN = x
+
     @SNOWFALL.setter
     def SNOWFALL(self, x):
         self.__SNOWFALL = x
+
     @LWin.setter
     def LWin(self, x):
         self.__LWin = x
+
     @LWout.setter
     def LWout(self, x):
         self.__LWout = x
+
     @H.setter
     def H(self, x):
         self.__H = x
+
     @LE.setter
     def LE(self, x):
         self.__LE = x
+
     @B.setter
     def B(self, x):
         self.__B = x
+
     @QRR.setter
     def QRR(self, x):
         self.__QRR = x
+
     @MB.setter
     def MB(self, x):
         self.__MB = x
@@ -968,7 +991,7 @@ def add_variable_along_latlon(
         Returns:
             Target dataset with the new spatial variable.
         """
-        ds[name] = ((Config.northing, Config.easting), var.data)
+        ds[name] = ((main_config.northing, main_config.easting), var.data)
         self.set_variable_metadata(ds[name], units, long_name)
         return ds
 
@@ -1006,7 +1029,7 @@ def add_variable_along_latlontime(
         Returns:
             Target dataset with the new spatiotemporal variable.
         """
-        ds[name] = (("time", Config.northing, Config.easting), var.data)
+        ds[name] = (("time", main_config.northing, main_config.easting), var.data)
         self.set_variable_metadata(ds[name], units, long_name)
         return ds
 
@@ -1025,7 +1048,7 @@ def add_variable_along_latlonlayertime(
         Returns:
             Target dataset with the new spatiotemporal mesh.
         """
-        ds[name] = (("time", Config.northing, Config.easting, "layer"), var.data)
+        ds[name] = (("time", main_config.northing, main_config.easting, "layer"), var.data)
         self.set_variable_metadata(ds[name], units, long_name)
         return ds
 
@@ -1044,7 +1067,7 @@ def add_variable_along_latlonlayer(
         Returns:
             Target dataset with the new spatial mesh.
         """
-        ds[name] = ((Config.northing, Config.easting, "layer"), var.data)
+        ds[name] = ((main_config.northing, main_config.easting, "layer"), var.data)
         self.set_variable_metadata(ds[name], units, long_name)
         return ds
 
diff --git a/cosipy/cpkernel/node.py b/cosipy/cpkernel/node.py
index 5527d17..f97ea0f 100644
--- a/cosipy/cpkernel/node.py
+++ b/cosipy/cpkernel/node.py
@@ -4,20 +4,20 @@
 from numba import float64
 from numba.experimental import jitclass
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config
 
 # only required for jitclass/njit
-ice_density = Constants.ice_density
-water_density = Constants.water_density
-air_density = Constants.air_density
-spec_heat_ice = Constants.spec_heat_ice
-spec_heat_air = Constants.spec_heat_air
-spec_heat_water = Constants.spec_heat_water
-k_i = Constants.k_i
-k_a = Constants.k_a
-k_w = Constants.k_w
-thermal_conductivity_method = Constants.thermal_conductivity_method
-zero_temperature = Constants.zero_temperature
+ice_density = constants_config.ice_density
+water_density = constants_config.water_density
+air_density = constants_config.air_density
+spec_heat_ice = constants_config.spec_heat_ice
+spec_heat_air = constants_config.spec_heat_air
+spec_heat_water = constants_config.spec_heat_water
+k_i = constants_config.k_i
+k_a = constants_config.k_a
+k_w = constants_config.k_w
+thermal_conductivity_method = constants_config.thermal_conductivity_method
+zero_temperature = constants_config.zero_temperature
 
 spec = OrderedDict()
 spec["height"] = float64
diff --git a/cosipy/modules/albedo.py b/cosipy/modules/albedo.py
index 2e6afd3..85b3d4f 100644
--- a/cosipy/modules/albedo.py
+++ b/cosipy/modules/albedo.py
@@ -1,21 +1,21 @@
 import numpy as np
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 # only required for njitted functions
-albedo_method = Constants.albedo_method
-albedo_mod_snow_aging = Constants.albedo_mod_snow_aging
-snow_ice_threshold = Constants.snow_ice_threshold
-albedo_firn = Constants.albedo_firn
-albedo_fresh_snow = Constants.albedo_fresh_snow
-albedo_ice = Constants.albedo_ice
-albedo_mod_snow_depth = Constants.albedo_mod_snow_depth
-dt = Constants.dt
-zero_temperature = Constants.zero_temperature
-t_star_cutoff = Constants.t_star_cutoff
-t_star_dry = Constants.t_star_dry
-t_star_wet = Constants.t_star_wet
-t_star_K = Constants.t_star_K
+albedo_method = cc.albedo_method
+albedo_mod_snow_aging = cc.albedo_mod_snow_aging
+snow_ice_threshold = cc.snow_ice_threshold
+albedo_firn = cc.albedo_firn
+albedo_fresh_snow = cc.albedo_fresh_snow
+albedo_ice = cc.albedo_ice
+albedo_mod_snow_depth = cc.albedo_mod_snow_depth
+dt = cc.dt
+zero_temperature = cc.zero_temperature
+t_star_cutoff = cc.t_star_cutoff
+t_star_dry = cc.t_star_dry
+t_star_wet = cc.t_star_wet
+t_star_K = cc.t_star_K
 
 
 def updateAlbedo(
diff --git a/cosipy/modules/densification.py b/cosipy/modules/densification.py
index 0189c0f..2610cf1 100644
--- a/cosipy/modules/densification.py
+++ b/cosipy/modules/densification.py
@@ -1,13 +1,13 @@
 import numpy as np
 from numba import njit
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 # only required for njitted functions
-densification_method = Constants.densification_method
-snow_ice_threshold = Constants.snow_ice_threshold
-minimum_snow_layer_height = Constants.minimum_snow_layer_height
-zero_temperature = Constants.zero_temperature
+densification_method = cc.densification_method
+snow_ice_threshold = cc.snow_ice_threshold
+minimum_snow_layer_height = cc.minimum_snow_layer_height
+zero_temperature = cc.zero_temperature
 
 
 def densification(GRID, SLOPE, dt):
@@ -259,7 +259,7 @@ def method_empirical(GRID, SLOPE, dt):
 
         if (1 - (dRho / GRID.get_node_density(idxNode))) < 1:
             # Set the new ice fraction
-            GRID.set_node_ice_fraction(idxNode, (rho_max + (rho[idxNode]-rho_max) * np.exp(-dt/tau))/Constants.ice_density )
+            GRID.set_node_ice_fraction(idxNode, (rho_max + (rho[idxNode]-rho_max) * np.exp(-dt/tau))/cc.ice_density )
 
             # Set height change
             GRID.set_node_height(idxNode, (1-(dRho/GRID.get_node_density(idxNode)))*GRID.get_node_height(idxNode))
diff --git a/cosipy/modules/evaluation.py b/cosipy/modules/evaluation.py
index e931190..3c4faed 100644
--- a/cosipy/modules/evaluation.py
+++ b/cosipy/modules/evaluation.py
@@ -1,4 +1,4 @@
-from cosipy.config import Config
+from cosipy.config import main_config
 
 
 def evaluate(stake_names, stake_data, df):
@@ -17,7 +17,7 @@ def evaluate(stake_names, stake_data, df):
         float or None: Statistical evaluation.
     """
 
-    if Config.eval_method == "rmse":
+    if main_config.eval_method == "rmse":
         stat = rmse(stake_names, stake_data, df)
     else:
         stat = None
@@ -36,12 +36,12 @@ def rmse(stake_names: list, stake_data, df) -> float:
     Returns:
         RMSE of simulated measurements.
     """
-    if Config.obs_type not in ["mb", "snowheight"]:
-        msg = f'RMSE not implemented for obs_type="{Config.obs_type}" in config.toml.'
+    if main_config.obs_type not in ["mb", "snowheight"]:
+        msg = f'RMSE not implemented for obs_type="{main_config.obs_type}" in config.toml.'
         raise NotImplementedError(msg)
     else:
         rmse = (
-            (stake_data[stake_names].subtract(df[Config.obs_type], axis=0))
+            (stake_data[stake_names].subtract(df[main_config.obs_type], axis=0))
             ** 2
         ).mean() ** 0.5
 
diff --git a/cosipy/modules/penetratingRadiation.py b/cosipy/modules/penetratingRadiation.py
index 0ce21f2..7763d6c 100644
--- a/cosipy/modules/penetratingRadiation.py
+++ b/cosipy/modules/penetratingRadiation.py
@@ -1,15 +1,15 @@
 import numpy as np
 from numba import njit
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 # only required for njitted functions
-zero_temperature = Constants.zero_temperature
-spec_heat_ice = Constants.spec_heat_ice
-ice_density = Constants.ice_density
-water_density = Constants.water_density
-snow_ice_threshold = Constants.snow_ice_threshold
-lat_heat_melting = Constants.lat_heat_melting
+zero_temperature = cc.zero_temperature
+spec_heat_ice = cc.spec_heat_ice
+ice_density = cc.ice_density
+water_density = cc.water_density
+snow_ice_threshold = cc.snow_ice_threshold
+lat_heat_melting = cc.lat_heat_melting
 
 
 def penetrating_radiation(GRID, SWnet: float, dt: int):
@@ -29,10 +29,10 @@ def penetrating_radiation(GRID, SWnet: float, dt: int):
         radiation.
     """
     penetrating_allowed = ["Bintanja95"]
-    if Constants.penetrating_method == "Bintanja95":
+    if cc.penetrating_method == "Bintanja95":
         subsurface_melt, Si = method_Bintanja(GRID, SWnet, dt)
     else:
-        msg = f'Penetrating method = "{Constants.penetrating_method}" is not allowed, must be one of {". ".join(penetrating_allowed)}'
+        msg = f'Penetrating method = "{cc.penetrating_method}" is not allowed, must be one of {". ".join(penetrating_allowed)}'
         raise ValueError(msg)
 
     return subsurface_melt, Si
diff --git a/cosipy/modules/refreezing.py b/cosipy/modules/refreezing.py
index b120374..d1a4624 100644
--- a/cosipy/modules/refreezing.py
+++ b/cosipy/modules/refreezing.py
@@ -1,15 +1,15 @@
 from numba import njit
 
-from cosipy.constants import Constants
-
-zero_temperature = Constants.zero_temperature
-air_density = Constants.air_density
-ice_density = Constants.ice_density
-water_density = Constants.water_density
-spec_heat_ice = Constants.spec_heat_ice
-spec_heat_water = Constants.spec_heat_water
-lat_heat_melting = Constants.lat_heat_melting
-snow_ice_threshold = Constants.snow_ice_threshold
+from cosipy.constants import constants_config as cc
+
+zero_temperature = cc.zero_temperature
+air_density = cc.air_density
+ice_density = cc.ice_density
+water_density = cc.water_density
+spec_heat_ice = cc.spec_heat_ice
+spec_heat_water = cc.spec_heat_water
+lat_heat_melting = cc.lat_heat_melting
+snow_ice_threshold = cc.snow_ice_threshold
 
 
 @njit
diff --git a/cosipy/modules/roughness.py b/cosipy/modules/roughness.py
index 4d90a97..6ccc1dd 100644
--- a/cosipy/modules/roughness.py
+++ b/cosipy/modules/roughness.py
@@ -1,4 +1,4 @@
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 
 def updateRoughness(GRID) -> float:
@@ -17,11 +17,11 @@ def updateRoughness(GRID) -> float:
     """
 
     roughness_allowed = ["Moelg12"]
-    if Constants.roughness_method == "Moelg12":
+    if cc.roughness_method == "Moelg12":
         sigma = method_Moelg(GRID)
     else:
         error_message = (
-            f'Roughness method = "{Constants.roughness_method}" is not allowed,',
+            f'Roughness method = "{cc.roughness_method}" is not allowed,',
             f'must be one of {", ".join(roughness_allowed)}',
         )
         raise ValueError(" ".join(error_message))
@@ -51,13 +51,13 @@ def method_Moelg(GRID) -> float:
     hours_since_snowfall = fresh_snow_timestamp / 3600.0
 
     # Check whether snow or ice
-    if GRID.get_node_density(0) <= Constants.snow_ice_threshold:
+    if GRID.get_node_density(0) <= cc.snow_ice_threshold:
         sigma = min(
-            Constants.roughness_fresh_snow
-            + Constants.aging_factor_roughness * hours_since_snowfall,
-            Constants.roughness_firn,
+            cc.roughness_fresh_snow
+            + cc.aging_factor_roughness * hours_since_snowfall,
+            cc.roughness_firn,
         )
     else:
-        sigma = Constants.roughness_ice  # Roughness length, set to ice
+        sigma = cc.roughness_ice  # Roughness length, set to ice
 
     return sigma / 1000
diff --git a/cosipy/modules/surfaceTemperature.py b/cosipy/modules/surfaceTemperature.py
index f952c0f..a79862d 100644
--- a/cosipy/modules/surfaceTemperature.py
+++ b/cosipy/modules/surfaceTemperature.py
@@ -2,24 +2,24 @@
 from numba import njit
 from scipy.optimize import minimize, newton
 
-from cosipy.config import Config
-from cosipy.constants import Constants
+from cosipy.config import main_config
+from cosipy.constants import constants_config as cc
 from cosipy.modules.secant import secant
 
-zlt1 = Constants.zlt1
-zlt2 = Constants.zlt2
-saturation_water_vapour_method = Constants.saturation_water_vapour_method
-zero_temperature = Constants.zero_temperature
-lat_heat_vaporize = Constants.lat_heat_vaporize
-lat_heat_sublimation = Constants.lat_heat_sublimation
-sigma = Constants.sigma
-stability_correction = Constants.stability_correction
-spec_heat_air = Constants.spec_heat_air
-spec_heat_water = Constants.spec_heat_water
-sfc_temperature_method = Constants.sfc_temperature_method
-surface_emission_coeff = Constants.surface_emission_coeff
-water_density = Constants.water_density
-WRF_X_CSPY = Config.WRF_X_CSPY
+zlt1 = cc.zlt1
+zlt2 = cc.zlt2
+saturation_water_vapour_method = cc.saturation_water_vapour_method
+zero_temperature = cc.zero_temperature
+lat_heat_vaporize = cc.lat_heat_vaporize
+lat_heat_sublimation = cc.lat_heat_sublimation
+sigma = cc.sigma
+stability_correction = cc.stability_correction
+spec_heat_air = cc.spec_heat_air
+spec_heat_water = cc.spec_heat_water
+sfc_temperature_method = cc.sfc_temperature_method
+surface_emission_coeff = cc.surface_emission_coeff
+water_density = cc.water_density
+WRF_X_CSPY = main_config.WRF_X_CSPY
 
 
 def update_surface_temperature(GRID, dt, z, z0, T2, rH2, p, SWnet, u2, RAIN, SLOPE, LWin=None, N=None):
diff --git a/cosipy/tests/test_GRID_getter_functions.py b/cosipy/tests/test_GRID_getter_functions.py
index b77e28c..a6a9010 100644
--- a/cosipy/tests/test_GRID_getter_functions.py
+++ b/cosipy/tests/test_GRID_getter_functions.py
@@ -2,7 +2,7 @@
 import pytest
 from numba import float64
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.grid import Grid
 
 
@@ -79,10 +79,10 @@ def get_ice_fraction(self, ice_fraction, density):
         if ice_fraction is None:
             a = (
                 density
-                - (1 - (density / Constants.ice_density))
-                * Constants.air_density
+                - (1 - (density / cc.ice_density))
+                * cc.air_density
             )
-            ice_fraction = a / Constants.ice_density
+            ice_fraction = a / cc.ice_density
         else:
             ice_fraction = ice_fraction
         return ice_fraction
@@ -93,10 +93,10 @@ def test_get_ice_fraction(self, arg_ice_fraction, conftest_boilerplate):
         if arg_ice_fraction is None:
             a = (
                 test_density
-                - (1 - (test_density / Constants.ice_density))
-                * Constants.air_density
+                - (1 - (test_density / cc.ice_density))
+                * cc.air_density
             )
-            test_ice = a / Constants.ice_density
+            test_ice = a / cc.ice_density
         else:
             test_ice = arg_ice_fraction
         compare_ice = self.get_ice_fraction(arg_ice_fraction, test_density)
@@ -188,7 +188,7 @@ def test_grid_get_number_snow_layers(self, grid, conftest_boilerplate):
             [
                 1
                 for idx in range(grid.number_nodes)
-                if grid.get_node_density(idx) < Constants.snow_ice_threshold
+                if grid.get_node_density(idx) < cc.snow_ice_threshold
             ]
         )
 
diff --git a/cosipy/tests/test_NODE_getter_functions.py b/cosipy/tests/test_NODE_getter_functions.py
index 9e3dc69..2039961 100644
--- a/cosipy/tests/test_NODE_getter_functions.py
+++ b/cosipy/tests/test_NODE_getter_functions.py
@@ -1,6 +1,6 @@
 import pytest
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.node import Node
 
 
@@ -107,9 +107,9 @@ def test_node_get_layer_ice_fraction(
         if arg_ice_fraction is None:
             test_ice_fraction = (
                 self.snow_density
-                - (1 - (self.snow_density / Constants.ice_density))
-                * Constants.air_density
-            ) / Constants.ice_density
+                - (1 - (self.snow_density / cc.ice_density))
+                * cc.air_density
+            ) / cc.ice_density
         else:
             test_ice_fraction = arg_ice_fraction
         compare_ice_fraction = node.get_layer_ice_fraction()
@@ -128,9 +128,9 @@ def test_node_get_layer_air_porosity(self, node, conftest_boilerplate):
 
     def test_node_get_layer_density(self, node, conftest_boilerplate):
         test_density = (
-            self.ice_fraction * Constants.ice_density
-            + self.lwc * Constants.water_density
-            + node.get_layer_air_porosity() * Constants.air_density
+            self.ice_fraction * cc.ice_density
+            + self.lwc * cc.water_density
+            + node.get_layer_air_porosity() * cc.air_density
         )
         assert conftest_boilerplate.check_output(
             node.get_layer_density(), float, test_density
@@ -144,9 +144,9 @@ def test_node_get_layer_porosity(self, node, conftest_boilerplate):
 
     def test_node_get_layer_specific_heat(self, node, conftest_boilerplate):
         test_specific_heat = (
-            (1 - self.lwc - self.ice_fraction) * Constants.spec_heat_air
-            + self.ice_fraction * Constants.spec_heat_ice
-            + self.lwc * Constants.spec_heat_water
+            (1 - self.lwc - self.ice_fraction) * cc.spec_heat_air
+            + self.ice_fraction * cc.spec_heat_ice
+            + self.lwc * cc.spec_heat_water
         )
         conftest_boilerplate.check_output(
             node.get_layer_specific_heat(), float, test_specific_heat
@@ -157,7 +157,7 @@ def test_node_get_layer_cold_content(self, node, conftest_boilerplate):
             -node.get_layer_specific_heat()
             * node.get_layer_density()
             * self.height
-            * (self.temperature - Constants.zero_temperature)
+            * (self.temperature - cc.zero_temperature)
         )
         conftest_boilerplate.check_output(
             node.get_layer_cold_content(), float, test_cold_content
@@ -167,9 +167,9 @@ def test_node_get_layer_thermal_conductivity(
         self, node, conftest_boilerplate
     ):
         test_thermal_conductivity = (
-            self.ice_fraction * Constants.k_i
-            + node.get_layer_porosity() * Constants.k_a
-            + self.lwc * Constants.k_w
+            self.ice_fraction * cc.k_i
+            + node.get_layer_porosity() * cc.k_a
+            + self.lwc * cc.k_w
         )
         conftest_boilerplate.check_output(
             node.get_layer_thermal_conductivity(),
diff --git a/cosipy/tests/test_config.py b/cosipy/tests/test_config.py
index 4e6589f..7775f65 100644
--- a/cosipy/tests/test_config.py
+++ b/cosipy/tests/test_config.py
@@ -1,14 +1,20 @@
 import argparse
+from pathlib import Path
+
+import pytest
 
 import cosipy.config
-from cosipy.config import Config
+from cosipy.config import Config, CosipyConfigModel, main_config
+
+file_path = Path("config.toml")
 
+@pytest.fixture()
+def test_cfg():
+    return Config(file_path).validate()
 
 class TestConfigParser:
     """Test argument parsing for config hook."""
 
-    file_path = "./config.toml"
-
     def test_set_parser(self):
         compare_parser = cosipy.config.set_parser()
         assert isinstance(compare_parser, argparse.ArgumentParser)
@@ -29,11 +35,9 @@ def test_set_parser(self):
 class TestConfig:
     """Test rtoml support."""
 
-    file_path = "config.toml"
-
-    def test_init_config(self):
-        test_cfg = Config()
-        assert isinstance(test_cfg, Config)
+    def test_init_config(self, test_cfg):
+        assert isinstance(test_cfg, CosipyConfigModel)
+        assert test_cfg == main_config
 
     def test_get_raw_toml(self):
         variable_names = [
@@ -49,16 +53,12 @@ def test_get_raw_toml(self):
             "SUBSET",
         ]
 
-        test_cfg = Config()
-        compare_toml = test_cfg.get_raw_toml()
+        compare_toml = Config(file_path).get_raw_toml(file_path)
         assert isinstance(compare_toml, dict)
         for name in variable_names:
-            assert name in compare_toml.keys()
-
-    def test_load_config(self):
-        test_cfg = Config()
-        test_cfg.load()
+            assert name in compare_toml
 
+    def test_load_config(self, test_cfg):
         variable_names = [
             "time_start",
             "time_end",
@@ -90,8 +90,3 @@ def test_load_config(self):
 
         for name in variable_names:
             assert hasattr(test_cfg, name)
-            compare_attr = getattr(test_cfg, name)
-            assert (
-                isinstance(compare_attr, (int, str, bool))
-                or compare_attr is None
-            )
diff --git a/cosipy/tests/test_parameterisation_albedo.py b/cosipy/tests/test_parameterisation_albedo.py
index 8975576..933d6f8 100644
--- a/cosipy/tests/test_parameterisation_albedo.py
+++ b/cosipy/tests/test_parameterisation_albedo.py
@@ -1,10 +1,11 @@
+import re
+
 import numpy as np
 import pytest
-import re
 
 import cosipy.modules.albedo as module_albedo
 from COSIPY import start_logging
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 
 class TestParamAlbedoUpdate:
@@ -16,13 +17,13 @@ def test_updateAlbedo(self, conftest_mock_grid):
         surface_albedo, snow_albedo = module_albedo.updateAlbedo(
             GRID=grid,
             surface_temperature=270.0,
-            albedo_snow=Constants.albedo_fresh_snow,
+            albedo_snow=cc.albedo_fresh_snow,
         )
         assert isinstance(surface_albedo, float)
         assert (
-            Constants.albedo_firn
+            cc.albedo_firn
             <= surface_albedo
-            <= Constants.albedo_fresh_snow
+            <= cc.albedo_fresh_snow
         )
         assert isinstance(snow_albedo, float)
 
@@ -33,22 +34,22 @@ def test_updateAlbedo_ice(
         albedo, snow_albedo = module_albedo.updateAlbedo(
             GRID=grid_ice,
             surface_temperature=270.0,
-            albedo_snow=Constants.albedo_fresh_snow,
+            albedo_snow=cc.albedo_fresh_snow,
         )
         assert conftest_boilerplate.check_output(
-            albedo, float, Constants.albedo_ice
+            albedo, float, cc.albedo_ice
         )
         assert conftest_boilerplate.check_output(
-            snow_albedo, float, Constants.albedo_fresh_snow
+            snow_albedo, float, cc.albedo_fresh_snow
         )
 
     @pytest.mark.parametrize("arg_height", [0.0, 0.05, 0.1, 0.5])
     @pytest.mark.parametrize(
         "arg_temperature",
-        [Constants.temperature_bottom, Constants.zero_temperature, 280.0],
+        [cc.temperature_bottom, cc.zero_temperature, 280.0],
     )
     @pytest.mark.parametrize(
-        "arg_time", [0, Constants.albedo_mod_snow_aging * 24]
+        "arg_time", [0, cc.albedo_mod_snow_aging * 24]
     )
     def test_get_surface_properties(
         self,
@@ -71,7 +72,7 @@ def test_get_surface_properties(
 
         grid.add_fresh_snow(
             arg_height,
-            Constants.constant_density,
+            cc.constant_density,
             arg_temperature,
             0.0,
             0,
@@ -105,7 +106,7 @@ def test_updateAlbedo_method(
         surface_albedo, snow_albedo = module_albedo.updateAlbedo(
             GRID=grid,
             surface_temperature=270.0,
-            albedo_snow=Constants.albedo_fresh_snow,
+            albedo_snow=cc.albedo_fresh_snow,
         )
         assert isinstance(surface_albedo, float)
         assert isinstance(snow_albedo, float)
@@ -133,7 +134,7 @@ def test_updateAlbedo_method_error(
             module_albedo.updateAlbedo(
                 GRID=grid,
                 surface_temperature=270.0,
-                albedo_snow=Constants.albedo_fresh_snow,
+                albedo_snow=cc.albedo_fresh_snow,
             )
 
 
@@ -144,9 +145,9 @@ def test_method_Oerlemans(self, conftest_mock_grid):
         """Get surface albedo without accounting for snow depth."""
 
         grid = conftest_mock_grid
-        albedo_limit = Constants.albedo_firn + (
-            Constants.albedo_fresh_snow - Constants.albedo_firn
-        ) * np.exp(0 / (Constants.albedo_mod_snow_aging * 24.0))
+        albedo_limit = cc.albedo_firn + (
+            cc.albedo_fresh_snow - cc.albedo_firn
+        ) * np.exp(0 / (cc.albedo_mod_snow_aging * 24.0))
 
         compare_albedo = module_albedo.method_Oerlemans(grid)
 
@@ -158,9 +159,9 @@ def test_method_Oerlemans(self, conftest_mock_grid):
     def test_get_simple_albedo(self, arg_hours):
         """Get surface albedo without accounting for snow depth."""
 
-        albedo_limit = Constants.albedo_firn + (
-            Constants.albedo_fresh_snow - Constants.albedo_firn
-        ) * np.exp((0) / (Constants.albedo_mod_snow_aging * 24.0))
+        albedo_limit = cc.albedo_firn + (
+            cc.albedo_fresh_snow - cc.albedo_firn
+        ) * np.exp((0) / (cc.albedo_mod_snow_aging * 24.0))
 
         compare_albedo = module_albedo.get_simple_albedo(
             elapsed_time=arg_hours
@@ -174,18 +175,18 @@ def test_get_simple_albedo(self, arg_hours):
     def test_get_albedo_with_decay(self, arg_depth):
         """Apply surface albedo decay due to the snow depth."""
 
-        albedo_limit = Constants.albedo_firn + (
-            Constants.albedo_fresh_snow - Constants.albedo_firn
-        ) * np.exp((0) / (Constants.albedo_mod_snow_aging * 24.0))
+        albedo_limit = cc.albedo_firn + (
+            cc.albedo_fresh_snow - cc.albedo_firn
+        ) * np.exp((0) / (cc.albedo_mod_snow_aging * 24.0))
 
-        albedo_limit = Constants.albedo_firn + (
-            Constants.albedo_ice - Constants.albedo_firn
+        albedo_limit = cc.albedo_firn + (
+            cc.albedo_ice - cc.albedo_firn
         ) * np.exp(
-            (-1.0 * arg_depth) / (Constants.albedo_mod_snow_depth / 100.0)
+            (-1.0 * arg_depth) / (cc.albedo_mod_snow_depth / 100.0)
         )
 
         compare_albedo = module_albedo.get_albedo_with_decay(
-            snow_albedo=Constants.albedo_firn, snow_height=arg_depth
+            snow_albedo=cc.albedo_firn, snow_height=arg_depth
         )
 
         assert isinstance(compare_albedo, float)
diff --git a/cosipy/tests/test_parameterisation_penetratingRadiation.py b/cosipy/tests/test_parameterisation_penetratingRadiation.py
index 469a087..bbbcecd 100644
--- a/cosipy/tests/test_parameterisation_penetratingRadiation.py
+++ b/cosipy/tests/test_parameterisation_penetratingRadiation.py
@@ -1,7 +1,7 @@
 import pytest
 
 import cosipy.modules.penetratingRadiation as pRad
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 
 class TestParamRadiation:
@@ -19,18 +19,18 @@ def test_penetrating_radiation_error(
         test_grid = conftest_mock_grid
         allow_list = ["Bintanja95"]
         conftest_boilerplate.patch_variable(
-            monkeypatch, pRad.Constants, {"penetrating_method": "WrongMethod"}
+            monkeypatch, pRad.cc, {"penetrating_method": "WrongMethod"}
         )
-        assert pRad.Constants.penetrating_method == "WrongMethod"
+        assert pRad.cc.penetrating_method == "WrongMethod"
         error_msg = (
-            f'Penetrating method = "{pRad.Constants.penetrating_method}" ',
+            f'Penetrating method = "{pRad.cc.penetrating_method}" ',
             f'is not allowed, must be one of {", ".join(allow_list)}',
         )
         with pytest.raises(ValueError, match="".join(error_msg)):
             pRad.penetrating_radiation(GRID=test_grid, SWnet=800.0, dt=3600)
 
     @pytest.mark.parametrize(
-        "arg_density", [250.0, Constants.snow_ice_threshold + 1]
+        "arg_density", [250.0, cc.snow_ice_threshold + 1]
     )
     def test_method_Bintanja(
         self,
@@ -40,19 +40,19 @@ def test_method_Bintanja(
         arg_density,
     ):
         conftest_boilerplate.patch_variable(
-            monkeypatch, pRad.Constants, {"penetrating_method": "Bintanja95"}
+            monkeypatch, pRad.cc, {"penetrating_method": "Bintanja95"}
         )
-        assert pRad.Constants.penetrating_method == "Bintanja95"
+        assert pRad.cc.penetrating_method == "Bintanja95"
         test_grid = conftest_mock_grid
         test_grid.add_fresh_snow(0.1, arg_density, 270.15, 0.0, 0.0)
         test_swnet = 800.0
-        if arg_density <= Constants.snow_ice_threshold:
+        if arg_density <= cc.snow_ice_threshold:
             test_si = test_swnet * 0.1
         else:
             test_si = test_swnet * 0.2
 
         melt_si = pRad.method_Bintanja(
-            GRID=test_grid, SWnet=test_swnet, dt=Constants.dt
+            GRID=test_grid, SWnet=test_swnet, dt=cc.dt
         )
         assert isinstance(melt_si, tuple)
         compare_melt = melt_si[0]
diff --git a/cosipy/tests/test_parameterisation_refreezing.py b/cosipy/tests/test_parameterisation_refreezing.py
index b85a5f9..6bb13e9 100644
--- a/cosipy/tests/test_parameterisation_refreezing.py
+++ b/cosipy/tests/test_parameterisation_refreezing.py
@@ -3,9 +3,9 @@
 import pytest
 
 from COSIPY import start_logging
+from cosipy.constants import constants_config as cc
 from cosipy.cpkernel.grid import Grid
 from cosipy.modules.refreezing import refreezing
-from cosipy.constants import Constants
 
 
 class TestParamRefreezing:
@@ -47,15 +47,15 @@ def test_percolation_convert_mwe(self, arg_mwe: float, arg_depth: float):
 
     def get_delta_theta_w(self, grid_obj):
         delta_temperature = (
-            np.array(grid_obj.get_temperature()) - Constants.zero_temperature
+            np.array(grid_obj.get_temperature()) - cc.zero_temperature
         )
 
         delta_theta_w = (
-            Constants.spec_heat_ice
-            * Constants.ice_density
+            cc.spec_heat_ice
+            * cc.ice_density
             * np.array(grid_obj.get_ice_fraction())
             * delta_temperature
-        ) / (Constants.lat_heat_melting * Constants.water_density)
+        ) / (cc.lat_heat_melting * cc.water_density)
 
         return delta_theta_w
 
@@ -69,7 +69,7 @@ def test_get_delta_theta_w(self, conftest_mock_grid):
     def get_delta_theta_i(self, theta_w):
 
         delta_theta_i = -(
-            Constants.water_density * theta_w / Constants.ice_density
+            cc.water_density * theta_w / cc.ice_density
         )
         return delta_theta_i
 
@@ -91,8 +91,8 @@ def test_refreezing(self):
         initial_ice_fraction = np.array(GRID.get_ice_fraction())
         initial_lwc = np.array(GRID.get_liquid_water_content())
         initial_mass = np.nansum(
-            initial_lwc * heights * Constants.water_density
-        ) + np.nansum(initial_ice_fraction * heights * Constants.ice_density)
+            initial_lwc * heights * cc.water_density
+        ) + np.nansum(initial_ice_fraction * heights * cc.ice_density)
 
         water_refreezed = refreezing(GRID)
         assert isinstance(water_refreezed, float)
@@ -104,21 +104,21 @@ def test_refreezing(self):
         final_mass = (
             np.nansum(
                 np.array(GRID.get_liquid_water_content())
-                * Constants.water_density
+                * cc.water_density
                 * heights
             )
             # + water_refreezed
             + np.nansum(
                 np.array(GRID.get_ice_fraction())
-                * Constants.ice_density
+                * cc.ice_density
                 * heights
             )
         )
 
         delta_w = final_lwc - initial_lwc
         delta_i = final_ice_fraction - initial_ice_fraction
-        delta_w_mass = -delta_w * Constants.water_density
-        delta_i_mass = delta_i * Constants.ice_density
+        delta_w_mass = -delta_w * cc.water_density
+        delta_i_mass = delta_i * cc.ice_density
 
         delta_i_mwe = delta_i_mass / heights
         delta_w_mwe = delta_w_mass / heights
diff --git a/cosipy/tests/test_parameterisation_roughness.py b/cosipy/tests/test_parameterisation_roughness.py
index 830db4c..662c728 100644
--- a/cosipy/tests/test_parameterisation_roughness.py
+++ b/cosipy/tests/test_parameterisation_roughness.py
@@ -1,9 +1,9 @@
 import numpy as np
 import pytest
 
-from COSIPY import start_logging
-from cosipy.constants import Constants
 import cosipy.modules.roughness as module_roughness
+from COSIPY import start_logging
+from cosipy.constants import constants_config as cc
 from cosipy.modules.roughness import updateRoughness
 
 
@@ -19,12 +19,12 @@ def test_updateRoughness_method_error(
 
         conftest_boilerplate.patch_variable(
             monkeypatch,
-            module_roughness.Constants,
+            module_roughness.cc,
             {"roughness_method": arg_method},
         )
-        assert module_roughness.Constants.roughness_method == arg_method
+        assert module_roughness.cc.roughness_method == arg_method
         error_message = (
-            f'Roughness method = "{module_roughness.Constants.roughness_method}"',
+            f'Roughness method = "{module_roughness.cc.roughness_method}"',
             f"is not allowed, must be one of",
             f'{", ".join(valid_methods)}',
         )
@@ -37,11 +37,11 @@ def test_method_moelg(self, conftest_mock_grid, conftest_mock_grid_ice):
         compare_roughness = module_roughness.method_Moelg(test_grid)
 
         assert (
-            Constants.roughness_fresh_snow / 1000
+            cc.roughness_fresh_snow / 1000
             <= compare_roughness
-            <= Constants.roughness_firn / 1000
+            <= cc.roughness_firn / 1000
         )
 
         test_grid_ice = conftest_mock_grid_ice
         ice_roughness = module_roughness.method_Moelg(test_grid_ice)
-        assert np.isclose(ice_roughness, Constants.roughness_ice / 1000)
+        assert np.isclose(ice_roughness, cc.roughness_ice / 1000)
diff --git a/cosipy/tests/test_parameterisation_surfaceTemperature.py b/cosipy/tests/test_parameterisation_surfaceTemperature.py
index 0eb9a51..8e0cbb9 100644
--- a/cosipy/tests/test_parameterisation_surfaceTemperature.py
+++ b/cosipy/tests/test_parameterisation_surfaceTemperature.py
@@ -2,7 +2,7 @@
 import pytest
 
 import cosipy.modules.surfaceTemperature as module_surface_temperature
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 
 
 class TestParamSurfaceTemperature:
@@ -41,7 +41,7 @@ class TestParamSurfaceTemperature:
     def assert_minimisation_ranges(
         self, surface_temperature, lw_in, lw_out, shf, lhf, ghf
     ):
-        assert Constants.zero_temperature >= surface_temperature >= 220.0
+        assert cc.zero_temperature >= surface_temperature >= 220.0
         assert 400 >= lw_in >= 0
         assert 0 >= lw_out >= -400
         assert 250 >= shf >= -250
@@ -56,12 +56,12 @@ def test_minimize_surface_energy_balance(
     ):
         conftest_boilerplate.patch_variable(
             monkeypatch,
-            Constants,
+            cc,
             {"sfc_temperature_method": arg_optim},
         )
         conftest_boilerplate.patch_variable(
             monkeypatch,
-            module_surface_temperature.Constants,
+            module_surface_temperature.cc,
             {"sfc_temperature_method": arg_optim},
         )
 
@@ -117,7 +117,7 @@ def test_minimize_newton(
     ):
         conftest_boilerplate.patch_variable(
             monkeypatch,
-            module_surface_temperature.Constants,
+            module_surface_temperature.cc,
             {"sfc_temperature_method": arg_optim},
         )
         test_grid = conftest_mock_grid
@@ -202,7 +202,7 @@ def test_surface_Temperature_parameterisation(self, conftest_mock_grid):
             N=0.5,  # otherwise tries to retrieve LWin from non-existent file
         )
 
-        assert Constants.zero_temperature >= surface_temperature >= 220.0
+        assert cc.zero_temperature >= surface_temperature >= 220.0
         assert 400 >= lw_radiation_in >= 0
         assert 0 >= lw_radiation_out >= -400
         assert 250 >= sensible_heat_flux >= -250
diff --git a/cosipy/tests/test_postprocessing_field_plots.py b/cosipy/tests/test_postprocessing_field_plots.py
index ed0d416..2ae1df9 100644
--- a/cosipy/tests/test_postprocessing_field_plots.py
+++ b/cosipy/tests/test_postprocessing_field_plots.py
@@ -5,6 +5,7 @@
 """
 
 import argparse
+from pathlib import Path
 from unittest.mock import patch
 
 import cartopy
@@ -361,7 +362,7 @@ def test_parse_arguments(self):
         assert isinstance(args, argparse.Namespace)
 
         # Required
-        assert args.file == "./path/file"
+        assert args.file == Path("./path/file")
         assert args.pdate == "'2009-01-01'"
         # Defaults
         assert args.variable is None
diff --git a/cosipy/utilities/aws2cosipy/aws2cosipy.py b/cosipy/utilities/aws2cosipy/aws2cosipy.py
index 6f875da..ef28c84 100644
--- a/cosipy/utilities/aws2cosipy/aws2cosipy.py
+++ b/cosipy/utilities/aws2cosipy/aws2cosipy.py
@@ -1,7 +1,6 @@
-"""
-Reads the input data (model forcing) and writes the output to a netCDF
-file. It supports point models with ``create_1D_input`` and distributed
-simulations with ``create_2D_input``.
+"""Reads the input data (model forcing) and writes the output to a netCDF file.
+It supports point models with ``create_1D_input`` and distributed simulations
+with ``create_2D_input``.
 
 The 1D input function works without a static file, in which the static
 variables are created.
@@ -67,8 +66,10 @@ def read_input_file(input_path: Path) -> pd.DataFrame:
     print(f"{'-' * 43}")
     print(f"Create input\nRead input file {input_path}")
 
-    date_parser = lambda x: dateutil.parser.parse(x, ignoretz=True)
-    dataframe = pd.read_csv(
+    def date_parser(x):
+        return dateutil.parser.parse(x, ignoretz=True)
+
+    return pd.read_csv(
         input_path,
         delimiter=_cfg.coords["delimiter"],
         index_col=["TIMESTAMP"],
@@ -77,13 +78,10 @@ def read_input_file(input_path: Path) -> pd.DataFrame:
         date_parser=date_parser,
     )
 
-    return dataframe
-
 
 def convert_to_numeric(series: pd.Series) -> pd.Series:
     """Convert series to numeric type."""
-    series = series.apply(pd.to_numeric, errors="coerce")
-    return series
+    return series.apply(pd.to_numeric, errors="coerce")
 
 
 def set_order_and_type(
@@ -115,10 +113,9 @@ def set_order_and_type(
     if (_cfg.names["LWin_var"] not in dataframe) and (
         _cfg.names["N_var"] not in dataframe
     ):
-        raise ValueError(
-            "No data for either incoming longwave radiation or cloud cover."
-        )
-    elif _cfg.names["LWin_var"] in dataframe:
+        msg = "No data for either incoming longwave radiation or cloud cover."
+        raise ValueError(msg)
+    if _cfg.names["LWin_var"] in dataframe:
         dataframe[_cfg.names["LWin_var"]] = convert_to_numeric(
             dataframe[_cfg.names["LWin_var"]]
         )
@@ -135,13 +132,13 @@ def set_order_and_type(
     return dataframe
 
 
-def check_data(dataset: xr.Dataset, dataframe: pd.DataFrame):
+def check_data(dataset: xr.Dataset, dataframe: pd.DataFrame) -> None:
     """Check data is within physically reasonable bounds."""
     check(dataset.T2, 316.16, 223.16)
     check(dataset.RH2, 100.0, 0.0)
     check(dataset.U2, 50.0, 0.0)
     check(dataset.G, 1600.0, 0.0)
-    check(dataset.PRES, 1080.0, 200.0)
+    check(dataset.PRESS, 1080.0, 200.0)
 
     if _cfg.names["RRR_var"] in dataframe:
         check(dataset.RRR, 20.0, 0.0)
@@ -153,7 +150,7 @@ def check_data(dataset: xr.Dataset, dataframe: pd.DataFrame):
         check(dataset.N, 1.0, 0.0)
 
 
-def get_time_slice(dataframe, start_date, end_date):
+def get_time_slice(dataframe: pd.DataFrame, start_date: str | None, end_date: str | None):
     if (start_date is not None) & (end_date is not None):
         dataframe = dataframe.loc[start_date:end_date]
     return dataframe
@@ -173,7 +170,6 @@ def set_bias(
         altitude: The height difference between a location and sensor.
         limit: If True, set negative values to zero. Default True.
     """
-
     biased_data = data + altitude * _cfg.lapse[lapse_type]
     if limit:
         biased_data = np.maximum(biased_data, 0.0)
@@ -183,7 +179,7 @@ def set_bias(
 
 def set_variable_metadata() -> dict:
     """Initialise variable names and units."""
-    metadata = {
+    return {
         "HGT": ("m", "Elevation"),
         "ASPECT": ("degrees", "Aspect of slope"),
         "SLOPE": ("degrees", "Terrain slope"),
@@ -192,14 +188,13 @@ def set_variable_metadata() -> dict:
         "RH2": ("%", "Relative humidity at 2 m"),
         "U2": ("m s\u207b\xb9", "Wind velocity at 2 m"),
         "G": ("W m\u207b\xb2", "Incoming shortwave radiation"),
-        "PRES": ("hPa", "Atmospheric Pressure"),
+        "PRESS": ("hPa", "Atmospheric Pressure"),
         "RRR": ("mm", "Total precipitation (liquid+solid)"),
         "SNOWFALL": ("m", "Snowfall"),
         "LWin": ("W m\u207b\xb2", "Incoming longwave radiation"),
         "N": ("%", "Cloud cover fraction"),
     }
 
-    return metadata
 
 
 def get_variable_metadata(name: str) -> dict:
@@ -215,7 +210,6 @@ def get_variable_metadata(name: str) -> dict:
 
 
 def set_dataset_coordinates(dataset, latitude, longitude):
-
     x, y = np.meshgrid(longitude, latitude)
     dataset.coords["lat"] = (("south_north", "west_east"), x)
     dataset.coords["lon"] = (("south_north", "west_east"), y)
@@ -223,9 +217,7 @@ def set_dataset_coordinates(dataset, latitude, longitude):
     return dataset
 
 
-def set_zero_field(
-    time_index: int, lat_index: int, lon_index: int
-) -> np.ndarray:
+def set_zero_field(time_index: int, lat_index: int, lon_index: int) -> np.ndarray:
     """Initialise array and fill with zeros."""
     field = np.zeros([time_index, lat_index, lon_index])
 
@@ -233,9 +225,7 @@ def set_zero_field(
 
 
 def get_pressure_bias(data, height):
-    slp = data / np.power(
-        (1 - (0.0065 * _cfg.station["stationAlt"]) / (288.15)), 5.255
-    )
+    slp = data / np.power((1 - (0.0065 * _cfg.station["stationAlt"]) / (288.15)), 5.255)
     pressure = slp * np.power((1 - (0.0065 * height) / (288.15)), 5.22)
 
     return pressure
@@ -264,7 +254,6 @@ def create_1D_input(
     Tobias Sauter 07.07.2019
     Anselm 04.07.2020
     """
-
     df = read_input_file(input_path=cs_file)
 
     print(df[pd.isnull(df).any(axis=1)])
@@ -376,7 +365,7 @@ def create_1D_input(
     U2 = df[_cfg.names["U2_var"]].values  # Wind velocity
     G = df[_cfg.names["G_var"]].values  # Incoming shortwave radiation
 
-    PRES = get_pressure_bias(
+    PRESS = get_pressure_bias(
         data=df[_cfg.names["PRES_var"]].values, height=_cfg.points["hgt"]
     )  # Pressure
 
@@ -421,7 +410,9 @@ def create_1D_input(
     RH2 = set_relative_humidity_bounds(RH2)
 
     # Add variables to file
-    add_variable_along_point(ds=ds, var=_cfg.points["hgt"], **get_variable_metadata("HGT"))
+    add_variable_along_point(
+        ds=ds, var=_cfg.points["hgt"], **get_variable_metadata("HGT")
+    )
     add_variable_along_point(ds=ds, var=aspect, **get_variable_metadata("ASPECT"))
     add_variable_along_point(ds=ds, var=slope, **get_variable_metadata("SLOPE"))
     add_variable_along_point(ds=ds, var=mask, **get_variable_metadata("MASK"))
@@ -429,14 +420,22 @@ def create_1D_input(
     add_variable_along_timelatlon_point(ds=ds, var=RH2, **get_variable_metadata("RH2"))
     add_variable_along_timelatlon_point(ds=ds, var=U2, **get_variable_metadata("U2"))
     add_variable_along_timelatlon_point(ds=ds, var=G, **get_variable_metadata("G"))
-    add_variable_along_timelatlon_point(ds=ds, var=PRES, **get_variable_metadata("PRES"))
+    add_variable_along_timelatlon_point(
+        ds=ds, var=PRESS, **get_variable_metadata("PRESS")
+    )
 
     if _cfg.names["RRR_var"] in df:
-        add_variable_along_timelatlon_point(ds=ds, var=RRR, **get_variable_metadata("RRR"))
+        add_variable_along_timelatlon_point(
+            ds=ds, var=RRR, **get_variable_metadata("RRR")
+        )
     if _cfg.names["SNOWFALL_var"] in df:
-        add_variable_along_timelatlon_point(ds, var=SNOWFALL, **get_variable_metadata("SNOWFALL"))
+        add_variable_along_timelatlon_point(
+            ds, var=SNOWFALL, **get_variable_metadata("SNOWFALL")
+        )
     if _cfg.names["LWin_var"] in df:
-        add_variable_along_timelatlon_point(ds=ds, var=LW, **get_variable_metadata("LWin"))
+        add_variable_along_timelatlon_point(
+            ds=ds, var=LW, **get_variable_metadata("LWin")
+        )
     if _cfg.names["N_var"] in df:
         add_variable_along_timelatlon_point(ds=ds, var=N, **get_variable_metadata("N"))
 
@@ -522,7 +521,7 @@ def create_2D_input(
     RH2 = df[_cfg.names["RH2_var"]]  # Relative humidity
     U2 = df[_cfg.names["U2_var"]]  # Wind velocity
     G = df[_cfg.names["G_var"]]  # Incoming shortwave radiation
-    PRES = df[_cfg.names["PRES_var"]]  # Pressure
+    PRESS = df[_cfg.names["PRES_var"]]  # Pressure
 
     if _cfg.names["in_K"]:
         T2 = df[_cfg.names["T2_var"]].values  # Temperature
@@ -569,13 +568,12 @@ def create_2D_input(
             data=RH2[t], lapse_type="lapse_RH", altitude=altitude, limit=False
         )
         U_interp[t, :, :] = U2[t]
+        """Interpolate pressure using the barometric equation.
 
-        """
-        Interpolate pressure using the barometric equation.
         Do not replace with get_pressure_bias() as the arrays won't
         interpolate correctly.
         """
-        slp = PRES[t] / np.power(
+        slp = PRESS[t] / np.power(
             (1 - (0.0065 * _cfg.station["stationAlt"]) / (288.15)), 5.255
         )
         P_interp[t, :, :] = slp * np.power(
@@ -694,9 +692,7 @@ def create_2D_input(
             print("Look-up-table for topographic shading done")
 
             # Save look-up tables
-            Nt = int(
-                366 * (3600 / _cfg.radiation["dtstep"]) * 24
-            )  # number of time steps
+            Nt = int(366 * (3600 / _cfg.radiation["dtstep"]) * 24)  # number of time steps
             Ny = len(lats)  # number of latitudes
             Nx = len(lons)  # number of longitudes
 
@@ -761,21 +757,31 @@ def create_2D_input(
 
     # Add variables to file
     add_variable_along_latlon(ds=dso, var=ds.HGT.values, **get_variable_metadata("HGT"))
-    add_variable_along_latlon(ds=dso, var=ds.ASPECT.values, **get_variable_metadata("ASPECT"))
-    add_variable_along_latlon(ds=dso, var=ds.SLOPE.values, **get_variable_metadata("SLOPE"))
+    add_variable_along_latlon(
+        ds=dso, var=ds.ASPECT.values, **get_variable_metadata("ASPECT")
+    )
+    add_variable_along_latlon(
+        ds=dso, var=ds.SLOPE.values, **get_variable_metadata("SLOPE")
+    )
     add_variable_along_latlon(ds=dso, var=ds.MASK.values, **get_variable_metadata("MASK"))
     add_variable_along_timelatlon(ds=dso, var=T_interp, **get_variable_metadata("T2"))
     add_variable_along_timelatlon(ds=dso, var=RH_interp, **get_variable_metadata("RH2"))
     add_variable_along_timelatlon(ds=dso, var=U_interp, **get_variable_metadata("U2"))
     add_variable_along_timelatlon(ds=dso, var=G_interp, **get_variable_metadata("G"))
-    add_variable_along_timelatlon(ds=dso, var=P_interp, **get_variable_metadata("PRES"))
+    add_variable_along_timelatlon(ds=dso, var=P_interp, **get_variable_metadata("PRESS"))
 
     if _cfg.names["RRR_var"] in df:
-        add_variable_along_timelatlon(ds=dso, var=RRR_interp, **get_variable_metadata("RRR"))
+        add_variable_along_timelatlon(
+            ds=dso, var=RRR_interp, **get_variable_metadata("RRR")
+        )
     if _cfg.names["SNOWFALL_var"] in df:
-        add_variable_along_timelatlon(ds=dso, var=SNOWFALL_interp, **get_variable_metadata("SNOWFALL"))
+        add_variable_along_timelatlon(
+            ds=dso, var=SNOWFALL_interp, **get_variable_metadata("SNOWFALL")
+        )
     if _cfg.names["LWin_var"] in df:
-        add_variable_along_timelatlon(ds=dso, var=LW_interp, **get_variable_metadata("LWin"))
+        add_variable_along_timelatlon(
+            ds=dso, var=LW_interp, **get_variable_metadata("LWin")
+        )
     if _cfg.names["N_var"] in df:
         add_variable_along_timelatlon(ds=dso, var=N_interp, **get_variable_metadata("N"))
 
@@ -839,7 +845,6 @@ def add_variable_along_point(ds, var, name, units, long_name):
 
 def check(field, max_bound, min_bound):
     """Check the validity of the input data."""
-
     if np.nanmax(field) > max_bound or np.nanmin(field) < min_bound:
         msg = f"{str.capitalize(field.name)} MAX: {np.nanmax(field):.2f} MIN: {np.nanmin(field):.2f}"
         print(
@@ -849,14 +854,10 @@ def check(field, max_bound, min_bound):
 
 def check_for_nan(ds):
     if _cfg.coords["WRF"] is True:
-        for y, x in product(
-            range(ds.dims["south_north"]), range(ds.dims["west_east"])
-        ):
+        for y, x in product(range(ds.dims["south_north"]), range(ds.dims["west_east"])):
             mask = ds.MASK.sel(south_north=y, west_east=x)
             if mask == 1:
-                if np.isnan(
-                    ds.sel(south_north=y, west_east=x).to_array()
-                ).any():
+                if np.isnan(ds.sel(south_north=y, west_east=x).to_array()).any():
                     raise_nan_error()
     else:
         for y, x in product(range(ds.dims["lat"]), range(ds.dims["lon"])):
@@ -880,15 +881,11 @@ def check_temperature_bounds(temperature: np.ndarray):
     min_temperature = np.nanmin(temperature)
     check_msg = "Please check the input temperature"
     if max_temperature > 373.16:
-        error_message = (
-            f"Maximum temperature is: {max_temperature} K. {check_msg}"
-        )
-        raise ValueError(error_message)
-    elif min_temperature < 173.16:
-        error_message = (
-            f"Minimum temperature is: {min_temperature} K. {check_msg}"
-        )
-        raise ValueError(error_message)
+        msg = f"Maximum temperature is: {max_temperature} K. {check_msg}"
+        raise ValueError(msg)
+    if min_temperature < 173.16:
+        msg = f"Minimum temperature is: {min_temperature} K. {check_msg}"
+        raise ValueError(msg)
 
 
 def set_relative_humidity_bounds(humidity):
@@ -916,8 +913,7 @@ def nansumwrapper(a, **kw_args):
     """Sum dataframe columns which contain NaNs."""
     if np.isnan(a).all():
         return np.nan
-    else:
-        return np.nansum(a, **kw_args)
+    return np.nansum(a, **kw_args)
 
 
 def get_user_arguments(parser: argparse.ArgumentParser) -> argparse.Namespace:
diff --git a/cosipy/utilities/config_utils.py b/cosipy/utilities/config_utils.py
index 090c0b2..17d4b3a 100644
--- a/cosipy/utilities/config_utils.py
+++ b/cosipy/utilities/config_utils.py
@@ -5,6 +5,9 @@
 import argparse
 import sys
 from collections import namedtuple
+from pathlib import Path
+
+from cosipy.config import default_utilities_path
 
 if sys.version_info >= (3, 11):
     import tomllib
@@ -16,7 +19,7 @@ class TomlLoader(object):
     """Load and parse configuration files."""
 
     @staticmethod
-    def get_raw_toml(file_path: str = "./utilities_config.toml") -> dict:
+    def get_raw_toml(file_path: Path = default_utilities_path) -> dict:
         """Open and load .toml configuration file.
 
         Args:
@@ -25,7 +28,7 @@ def get_raw_toml(file_path: str = "./utilities_config.toml") -> dict:
         Returns:
             Loaded .toml data.
         """
-        with open(file_path, "rb") as f:
+        with file_path.open("rb") as f:
             raw_config = tomllib.load(f)
 
         return raw_config
@@ -43,11 +46,11 @@ class UtilitiesConfig(TomlLoader):
         wrf2cosipy: Configuration parameters for `wrf2cosipy.py`.
     """
 
-    def __init__(self):
+    def __init__(self) -> None:
         self.parser = self.set_arg_parser()
 
     @classmethod
-    def load(cls, path: str = "./utilities_config.toml"):
+    def load(cls, path: Path = default_utilities_path) -> None:
         raw_toml = cls.get_raw_toml(path)
         cls.set_config_values(raw_toml)
 
@@ -75,7 +78,7 @@ def set_arg_parser(cls) -> argparse.ArgumentParser:
             "--utilities",
             default="./utilities_config.toml",
             dest="utilities_path",
-            type=str,
+            type=Path,
             metavar="<path>",
             required=False,
             help="relative path to utilities' configuration file",
diff --git a/cosipy/utilities/wrf2cosipy/wrf2cosipy.py b/cosipy/utilities/wrf2cosipy/wrf2cosipy.py
index 8a1efd1..0200e16 100644
--- a/cosipy/utilities/wrf2cosipy/wrf2cosipy.py
+++ b/cosipy/utilities/wrf2cosipy/wrf2cosipy.py
@@ -34,7 +34,7 @@
 import xarray as xr
 import xrspatial as xrs
 
-from cosipy.constants import Constants
+from cosipy.constants import constants_config as cc
 from cosipy.utilities.config_utils import UtilitiesConfig
 
 _args = None
@@ -99,16 +99,16 @@ def create_input(wrf_file: Path, cosipy_file: Path, start_date, end_date):
     dso = add_variable_along_timelatlon(dso, rrr, 'RRR', 'mm', 'Total precipitation')
  
     # Calc fresh snow density
-    if (Constants.densification_method!='constant'):
+    if (cc.densification_method!='constant'):
         density_fresh_snow = np.maximum(109.0+6.0*(ds.T2.values-273.16)+26.0*np.sqrt(U2), 50.0)
     else:
-        density_fresh_snow = Constants.constant_density
+        density_fresh_snow = cc.constant_density
 
     # Derive snowfall from accumulated values
     snowf = np.full_like(ds.T2, 0.)
     for t in np.arange(1,len(dso.time)):
         snowf[t,:,:] = ds.SNOWNC[t,:,:]-ds.SNOWNC[t-1,:,:]
-    snowf = (snowf/1000.0)*(Constants.water_density/density_fresh_snow)
+    snowf = (snowf/1000.0)*(cc.water_density/density_fresh_snow)
     dso = add_variable_along_timelatlon(dso, snowf, 'SNOWFALL', 'm', 'Snowfall')
 
     # Compute slope & aspect from HGT
diff --git a/dev_requirements.txt b/dev_requirements.txt
index 0f02525..b7bf689 100644
--- a/dev_requirements.txt
+++ b/dev_requirements.txt
@@ -16,6 +16,7 @@ nco
 cdo
 cartopy
 vtk
+pydantic
 richdem; python_version < '3.11'
 tomli; python_version < '3.11'
 coveralls
diff --git a/docs/requirements.txt b/docs/requirements.txt
index 3b4ab97..cbd26fa 100644
--- a/docs/requirements.txt
+++ b/docs/requirements.txt
@@ -14,6 +14,7 @@ nco
 cdo
 cartopy
 vtk
+pydantic
 richdem; python_version < '3.11'
 pytest
 pytest-dependency
diff --git a/requirements.txt b/requirements.txt
index 189d950..6aad901 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -15,6 +15,7 @@ nco
 cdo
 cartopy
 vtk
+pydantic
 richdem; python_version < '3.11'
 coveralls
 codecov