fix wrong graupel tuning parameter config in data test

model/atmosphere/dycore/tests/dycore/fixtures.py

@@ -18,6 +18,7 @@
     grid_savepoint,
     htop_moist_proc,
     icon_grid,
+    icon_namelist,
     interpolation_savepoint,
     istep_exit,
     istep_init,
@@ -38,9 +39,11 @@
     savepoint_velocity_exit,
     savepoint_velocity_init,
     savepoint_vertically_implicit_dycore_solver_init,
+    solve_nonhydro_config,
     step_date_exit,
     step_date_init,
     stretch_factor,
     substep_exit,
     substep_init,
+    vertical_grid_config,
 )

model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/microphysics_options.py

@@ -5,6 +5,8 @@
 #
 # Please, refer to the LICENSE file in the root directory.
 # SPDX-License-Identifier: BSD-3-Clause
+from typing import Literal
+
 import gt4py.next as gtx
 from gt4py.eve import utils as eve_utils
 
@@ -29,3 +31,12 @@ class SnowInterceptParametererization(eve_utils.FrozenNamespace[gtx.int32]):
     FIELD_BEST_FIT_ESTIMATION = 1
     #: Estimated intercept parameter for the snow size distribution from the general moment equation in table 2 of Field et al. (2005)
     FIELD_GENERAL_MOMENT_ESTIMATION = 2
+
+
+ValidLiquidAutoConversionType = Literal[
+    LiquidAutoConversionType.KESSLER, LiquidAutoConversionType.SEIFERT_BEHENG
+]
+ValidSnowInterceptParametererization = Literal[
+    SnowInterceptParametererization.FIELD_BEST_FIT_ESTIMATION,
+    SnowInterceptParametererization.FIELD_GENERAL_MOMENT_ESTIMATION,
+]

model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/single_moment_six_class_gscp_graupel.py

@@ -61,19 +61,19 @@ class SingleMomentSixClassIconGraupelConfig:
     """
 
     #: liquid auto conversion mode. Originally defined as isnow_n0temp (PARAMETER) in gscp_data.f90 in ICON. I keep it because I think the choice depends on resolution.
-    liquid_autoconversion_option: mphys_options.LiquidAutoConversionType = (
+    liquid_autoconversion_option: mphys_options.ValidLiquidAutoConversionType = (
         mphys_options.LiquidAutoConversionType.KESSLER
     )
     #: snow size distribution interception parameter. Originally defined as isnow_n0temp (PARAMETER) in gscp_data.f90 in ICON. I keep it because I think the choice depends on resolution.
-    snow_intercept_option: mphys_options.SnowInterceptParametererization = (
+    snow_intercept_option: mphys_options.ValidSnowInterceptParametererization = (
         mphys_options.SnowInterceptParametererization.FIELD_GENERAL_MOMENT_ESTIMATION
     )
     #: Do latent heat nudging. Originally defined as dass_lhn in mo_run_config.f90 in ICON.
     do_latent_heat_nudging = False
     #: Whether a fixed latent heat capacities are used for water. Originally defined as ithermo_water in mo_nwp_tuning_config.f90 in ICON (0 means True).
     use_constant_latent_heat = True
     #: First parameter in RHS of eq. 5.163 in the COSMO microphysics documentation for the sticking efficiency when lstickeff = True (repricated in icon4py because it is always True in ICON). Originally defined as tune_zceff_min in mo_tuning_nwp_config.f90 in ICON.
-    ice_stickeff_min: ta.wpfloat = 0.075
+    ice_stickeff_min: ta.wpfloat = 0.01
     #: Power law coefficient in v-qi ice terminal velocity-mixing ratio relationship, see eq. 5.169 in the COSMO microphysics documentation. Originally defined as tune_zvz0i in mo_tuning_nwp_config.f90 in ICON.
     power_law_coeff_for_ice_mean_fall_speed: ta.wpfloat = 1.25
     #: Exponent of the density factor in ice terminal velocity equation to account for density (air thermodynamic state) change. Originally defined as tune_icesedi_exp in mo_tuning_nwp_config.f90 in ICON.

model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/stencils/microphyiscal_processes.py

@@ -31,7 +31,7 @@
 
 @gtx.field_operator
 def compute_cooper_inp_concentration(temperature: ta.wpfloat) -> ta.wpfloat:
-    cnin = 5.0 * exp(0.304 * (_phy_const.tmelt - temperature))
+    cnin = wpfloat("5.0") * exp(wpfloat("0.304") * (_phy_const.tmelt - temperature))
     cnin = minimum(cnin, _microphy_const.NIMAX_THOM)
     return cnin
 

model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_saturation_adjustment.py

@@ -33,32 +33,29 @@
 @pytest.mark.embedded_static_args
 @pytest.mark.datatest
 @pytest.mark.parametrize(
-    "experiment, model_top_height",
-    [(definitions.Experiments.WEISMAN_KLEMP_TORUS, 30000.0)],
+    "experiment",
+    [definitions.Experiments.WEISMAN_KLEMP_TORUS],
 )
 @pytest.mark.parametrize(
     "date", ["2008-09-01T01:59:48.000", "2008-09-01T01:59:52.000", "2008-09-01T01:59:56.000"]
 )
 @pytest.mark.parametrize("location", ["nwp-gscp-interface", "interface-nwp"])
 def test_saturation_adjustement(
-    location: str,
-    model_top_height: ta.wpfloat,
     date: str,
+    location: str,
     *,
     data_provider: sb.IconSerialDataProvider,
     grid_savepoint: sb.IconGridSavepoint,
     metrics_savepoint: sb.MetricSavepoint,
     icon_grid: base_grid.Grid,
+    model_top_height: ta.wpfloat,
     backend: gtx_typing.Backend,
 ) -> None:
     satad_init = data_provider.from_savepoint_satad_init(location=location, date=date)
     satad_exit = data_provider.from_savepoint_satad_exit(location=location, date=date)
 
-    config = satad.SaturationAdjustmentConfig(
-        tolerance=1e-3,
-        max_iter=10,
-    )
-    dtime = 2.0
+    config = satad.SaturationAdjustmentConfig()
+    dtime = data_provider.from_savepoint_weisman_klemp_graupel_entry(date=date).dtime()
 
     vertical_config = v_grid.VerticalGridConfig(
         icon_grid.num_levels,

model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_single_moment_six_class_gscp_graupel.py

@@ -38,30 +38,27 @@
 @pytest.mark.embedded_static_args
 @pytest.mark.datatest
 @pytest.mark.parametrize(
-    "experiment, model_top_height",
+    "experiment",
     [
-        (definitions.Experiments.WEISMAN_KLEMP_TORUS, 30000.0),
+        definitions.Experiments.WEISMAN_KLEMP_TORUS,
     ],
 )
 @pytest.mark.parametrize(
     "date", ["2008-09-01T01:59:48.000", "2008-09-01T01:59:52.000", "2008-09-01T01:59:56.000"]
 )
 def test_graupel(
     experiment: definitions.Experiments,
-    model_top_height: ta.wpfloat,
     date: str,
     *,
     data_provider: sb.IconSerialDataProvider,
     grid_savepoint: sb.IconGridSavepoint,
     metrics_savepoint: sb.MetricSavepoint,
     icon_grid: icon_grid.IconGrid,
-    lowest_layer_thickness: ta.wpfloat,
+    model_top_height: ta.wpfloat,
     backend: gtx_typing.Backend,
 ):
-    pytest.xfail("Tolerances have increased with new ser_data, need to check with @ongchia")
     vertical_config = v_grid.VerticalGridConfig(
         icon_grid.num_levels,
-        lowest_layer_thickness=lowest_layer_thickness,
         model_top_height=model_top_height,
     )
     vertical_params = v_grid.VerticalGrid(
@@ -103,16 +100,7 @@ def test_graupel(
         v=None,
     )
 
-    graupel_config = graupel.SingleMomentSixClassIconGraupelConfig(
-        liquid_autoconversion_option=mphys_options.LiquidAutoConversionType.SEIFERT_BEHENG,
-        ice_stickeff_min=0.01,
-        power_law_coeff_for_ice_mean_fall_speed=1.25,
-        exponent_for_density_factor_in_ice_sedimentation=0.3,
-        power_law_coeff_for_snow_fall_speed=20.0,
-        rain_mu=0.0,
-        rain_n0=1.0,
-        snow2graupel_riming_coeff=0.5,
-    )
+    graupel_config = definitions.construct_gscp_graupel_config(experiment)
 
     graupel_microphysics = graupel.SingleMomentSixClassIconGraupel(
         graupel_config=graupel_config,

model/testing/pyproject.toml

@@ -33,7 +33,8 @@ dependencies = [
   "pytest>=8.0.1",
   "serialbox4py>=2.6.2",
   "typing-extensions>=4.11.0",
-  "wget>=3.2"
+  "wget>=3.2",
+  "f90nml>=1.5"
 ]
 description = "Testing utils for the icon4py model."
 license = {text = "BSD-3 License"}
@@ -71,6 +72,12 @@ version = "{current_version}"
 [[tool.bumpversion.files]]
 filename = "src/icon4py/model/testing/__init__.py"
 
+# -- mypy --
+[tool.mypy]
+[[tool.mypy.overrides]]
+ignore_missing_imports = true
+module = ["f90nml.*"]
+
 # -- ruff --
 [tool.ruff]
 extend = "../../pyproject.toml"

model/testing/src/icon4py/model/testing/datatest_utils.py

@@ -11,6 +11,7 @@
 import pathlib
 import urllib.parse
 
+import f90nml
 import gt4py.next.typing as gtx_typing
 
 from icon4py.model.common.decomposition import definitions as decomposition
@@ -72,3 +73,15 @@ def create_icon_serial_data_provider(
         mpi_rank=rank,
         do_print=True,
     )
+
+
+def read_namelist(namelist_path: pathlib.Path) -> dict:
+    """Read namelist file and return a dictionary of parameters."""
+    namelist = f90nml.read(namelist_path)
+    # TODO (Chia Rui): This needs to be checked in the review whether we can just keep the first element of lists for all namelist parameters.
+    for namelist_name, namelist_content in namelist.items():
+        if isinstance(namelist_content, dict):
+            for variable_name, variable_value in namelist_content.items():
+                if isinstance(variable_value, list):
+                    namelist[namelist_name][variable_name] = variable_value[0]
+    return namelist

model/testing/src/icon4py/model/testing/definitions.py

@@ -20,11 +20,15 @@
 if TYPE_CHECKING:
     from icon4py.model.atmosphere.diffusion import diffusion
     from icon4py.model.atmosphere.dycore import solve_nonhydro as solve_nh
+    from icon4py.model.atmosphere.subgrid_scale_physics.microphysics import (
+        single_moment_six_class_gscp_graupel as graupel,
+    )
 
 
 SERIALIZED_DATA_DIR: Final = "ser_icondata"
 SERIALIZED_DATA_SUBDIR: Final = "ser_data"
 GRID_DATA_DIR: Final = "grids"
+NAMELIST_FILENAME: Final = "NAMELIST_ICON_output_atm"
 
 
 def serialized_data_path() -> pathlib.Path:
@@ -198,6 +202,7 @@ def construct_diffusion_config(
 ) -> diffusion.DiffusionConfig:
     from icon4py.model.atmosphere.diffusion import diffusion
 
+
     if experiment == Experiments.MCH_CH_R04B09:
         return diffusion.DiffusionConfig(
             diffusion_type=diffusion.DiffusionType.SMAGORINSKY_4TH_ORDER,
@@ -264,6 +269,25 @@ def construct_nonhydrostatic_config(experiment: Experiment) -> solve_nh.NonHydro
         )
 
 
+def construct_gscp_graupel_config(
+    experiment: Experiment,
+) -> graupel.SingleMomentSixClassIconGraupelConfig:
+    from icon4py.model.atmosphere.subgrid_scale_physics.microphysics import (
+        microphysics_options as mphys_options,
+        single_moment_six_class_gscp_graupel as graupel,
+    )
+
+    if experiment == Experiments.WEISMAN_KLEMP_TORUS:
+        return graupel.SingleMomentSixClassIconGraupelConfig(
+            liquid_autoconversion_option=mphys_options.LiquidAutoConversionType.SEIFERT_BEHENG,
+            ice_stickeff_min=0.075,
+        )
+    else:
+        raise NotImplementedError(
+            f"SingleMomentSixClassIconGraupelConfig for experiment {experiment.name} not implemented."
+        )
+
+
 def construct_metrics_config(experiment: Experiment) -> tuple:
     match experiment:
         case Experiments.MCH_CH_R04B09:
@@ -301,13 +325,13 @@ def construct_metrics_config(experiment: Experiment) -> tuple:
             thhgtd_zdiffu = 200.0
         case Experiments.WEISMAN_KLEMP_TORUS:
             lowest_layer_thickness = 50.0
-            model_top_height = 23500.0
-            stretch_factor = 1.0
-            damping_height = 8000.0
+            model_top_height = 30000.0
+            stretch_factor = 0.85
+            damping_height = 12000.0
             rayleigh_coeff = 0.75
             exner_expol = 0.333
-            vwind_offctr = 0.15
-            rayleigh_type = 2
+            vwind_offctr = 0.20
+            rayleigh_type = 3
             thslp_zdiffu = 0.025
             thhgtd_zdiffu = 125.0
         case _:

model/testing/src/icon4py/model/testing/fixtures/datatest.py

@@ -14,9 +14,10 @@
 import pytest
 
 import icon4py.model.common.decomposition.definitions as decomposition
+from icon4py.model.atmosphere.dycore import solve_nonhydro as solve_nh
 from icon4py.model.common import model_backends, model_options
 from icon4py.model.common.constants import RayleighType
-from icon4py.model.common.grid import base as base_grid
+from icon4py.model.common.grid import base as base_grid, vertical as v_grid
 from icon4py.model.testing import data_handling, datatest_utils as dt_utils, definitions
 
 
@@ -154,6 +155,74 @@ def data_provider(
     return dt_utils.create_icon_serial_data_provider(data_path, processor_props.rank, backend)
 
 
+@pytest.fixture
+def icon_namelist(
+    download_ser_data: None,  # downloads data as side-effect
+    experiment: definitions.Experiment,
+    processor_props: decomposition.ProcessProperties,
+) -> dict:
+    experiment_dir = dt_utils.get_ranked_experiment_name_with_version(
+        experiment,
+        processor_props.comm_size,
+    )
+    namelist_path = definitions.serialized_data_path().joinpath(
+        experiment_dir, definitions.NAMELIST_FILENAME
+    )
+    return dt_utils.read_namelist(namelist_path)
+
+
+@pytest.fixture
+def vertical_grid_config(
+    icon_namelist: dict,
+) -> v_grid.VerticalGridConfig:
+    return v_grid.VerticalGridConfig(
+        # TODO (Chia Rui): where should we put the config construction? And remove this hardcoded style in next commits
+        num_levels=icon_namelist["RUN_NML"]["NUM_LEV"],
+        maximal_layer_thickness=icon_namelist["SLEVE_NML"]["MAX_LAY_THCKN"],
+        top_height_limit_for_maximal_layer_thickness=icon_namelist["SLEVE_NML"]["HTOP_THCKNLIMIT"],
+        lowest_layer_thickness=icon_namelist["SLEVE_NML"]["MIN_LAY_THCKN"],
+        model_top_height=icon_namelist["SLEVE_NML"]["TOP_HEIGHT"],
+        flat_height=icon_namelist["SLEVE_NML"]["FLAT_HEIGHT"],
+        stretch_factor=icon_namelist["SLEVE_NML"]["STRETCH_FAC"],
+        rayleigh_damping_height=icon_namelist["NONHYDROSTATIC_NML"]["DAMP_HEIGHT"],
+        htop_moist_proc=icon_namelist["NONHYDROSTATIC_NML"]["HTOP_MOIST_PROC"],
+        SLEVE_decay_scale_1=icon_namelist["SLEVE_NML"]["DECAY_SCALE_1"],
+        SLEVE_decay_scale_2=icon_namelist["SLEVE_NML"]["DECAY_SCALE_2"],
+        SLEVE_decay_exponent=icon_namelist["SLEVE_NML"]["DECAY_EXP"],
+    )
+
+
+@pytest.fixture
+def solve_nonhydro_config(
+    icon_namelist: dict,
+) -> solve_nh.NonHydrostaticConfig:
+    return solve_nh.NonHydrostaticConfig(
+        itime_scheme=icon_namelist["NONHYDROSTATIC_NML"]["ITIME_SCHEME"],
+        iadv_rhotheta=icon_namelist["NONHYDROSTATIC_NML"]["IADV_RHOTHETA"],
+        igradp_method=icon_namelist["NONHYDROSTATIC_NML"]["IGRADP_METHOD"],
+        rayleigh_type=icon_namelist["NONHYDROSTATIC_NML"]["RAYLEIGH_TYPE"],
+        rayleigh_coeff=icon_namelist["NONHYDROSTATIC_NML"]["RAYLEIGH_COEFF"],
+        divdamp_order=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_ORDER"],
+        is_iau_active=False,  # TODO (Chia RUi): a bug to be fixed in https://github.com/C2SM/icon4py/pull/972
+        iau_wgt_dyn=0.0,  # TODO (Chia RUi): a bug to be fixed in https://github.com/C2SM/icon4py/pull/972
+        divdamp_type=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_TYPE"],
+        divdamp_trans_start=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_TRANS_START"],
+        divdamp_trans_end=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_TRANS_END"],
+        l_vert_nested=icon_namelist["RUN_NML"]["LVERT_NEST"],
+        rhotheta_offctr=icon_namelist["NONHYDROSTATIC_NML"]["RHOTHETA_OFFCTR"],
+        veladv_offctr=icon_namelist["NONHYDROSTATIC_NML"]["VELADV_OFFCTR"],
+        max_nudging_coefficient=icon_namelist["INTERPOL_NML"]["NUDGE_MAX_COEFF"],
+        fourth_order_divdamp_factor=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_FAC"],
+        fourth_order_divdamp_factor2=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_FAC2"],
+        fourth_order_divdamp_factor3=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_FAC3"],
+        fourth_order_divdamp_factor4=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_FAC4"],
+        fourth_order_divdamp_z=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_Z"],
+        fourth_order_divdamp_z2=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_Z2"],
+        fourth_order_divdamp_z3=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_Z3"],
+        fourth_order_divdamp_z4=icon_namelist["NONHYDROSTATIC_NML"]["DIVDAMP_Z4"],
+    )
+
+
 @pytest.fixture
 def grid_savepoint(
     data_provider: serialbox.IconSerialDataProvider, experiment: definitions.Experiment