fix wrong graupel tuning parameter config in data test

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/src/icon4py/model/testing/datatest_utils.py

@@ -10,6 +10,9 @@
 
 import pathlib
 import urllib.parse
+import ast
+import re
+
 
 import gt4py.next.typing as gtx_typing
 
@@ -72,3 +75,85 @@ def create_icon_serial_data_provider(
         mpi_rank=rank,
         do_print=True,
     )
+
+
+def read_namelist(path: pathlib.Path) -> dict:
+    """ICON NAMELIST_ICON_output_atm reader.
+    Returns a dictionary of dictionaries, where the keys are the namelist names
+    and the values are dictionaries of key-value pairs from the namelist.
+
+    Use as:
+    namelists = read_namelist("/path/to/NAMELIST_ICON_output_atm")
+    print(namelists["NWP_TUNING_NML"]["TUNE_ZCEFF_MIN"])
+    """
+    with path.open() as f:
+        txt = f.read()
+    namelist_set = re.findall(r"&(\w+)(.*?)\/", txt, re.S)
+    full_namelist = {}
+    for namelist_name, namelist_content in namelist_set:
+        full_namelist[namelist_name] = _parse_namelist_content(namelist_content)
+    return full_namelist
+
+
+def _parse_namelist_content(namelist_content):
+    """
+    Parse the contents of a single namelist group to a Python dictionary.
+    """
+    result = {}
+    current_variable = None
+
+    # Remove comments
+    namelist_content = re.sub(r"!.*", "", namelist_content)
+
+    # Split into lines
+    lines = namelist_content.splitlines()
+
+    for line in lines:
+        line = line.strip()
+        # skip any non-meaningful empty line
+        if not line:
+            continue
+
+        # Remove trailing comma
+        if line.endswith(","):
+            line = line[:-1]
+
+        if "=" in line:
+            # New variable-value pair
+            variable, value = line.split("=", 1)
+            current_variable = variable.strip()
+            result[current_variable] = _parse_value(value)
+        # TODO(Chia Rui): check if continuation lines is irrelevant in our tests
+        # else:
+        #     # Continuation line (array or multi-line string)
+        #     if current_variable is not None:
+        #         val = _parse_value(line)
+        #         # convert to a list if we have multiple values for the same variable
+        #         if not isinstance(result[current_variable], list):
+        #             result[current_variable] = [result[current_variable]]
+        #         result[current_variable].append(val)
+
+    return result
+
+
+def _parse_value(value):
+    """
+    Convert Fortran-style values to Python types.
+    """
+    value = value.strip()
+
+    # Fortran logical
+    if value in ("T", ".T."):
+        return True
+    if value in ("F", ".F."):
+        return False
+
+    # Quoted string
+    if value.startswith("'") and value.endswith("'"):
+        return value[1:-1].rstrip()
+
+    # Try numeric conversion
+    try:
+        return ast.literal_eval(value)
+    except Exception:
+        return value

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

@@ -16,7 +16,7 @@
 import icon4py.model.common.decomposition.definitions as decomposition
 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 +154,42 @@ 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 grid_savepoint(
     data_provider: serialbox.IconSerialDataProvider, experiment: definitions.Experiment