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/definitions.py

@@ -20,6 +20,9 @@
 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"
@@ -264,6 +267,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 +323,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/namelist_reader.py

@@ -0,0 +1,39 @@
+# ICON4Py - ICON inspired code in Python and GT4Py
+#
+# Copyright (c) 2022-2024, ETH Zurich and MeteoSwiss
+# All rights reserved.
+#
+# Please, refer to the LICENSE file in the root directory.
+# SPDX-License-Identifier: BSD-3-Clause
+import ast
+import pathlib
+import re
+
+
+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()
+    blocks = re.findall(r"&(\w+)(.*?)\/", txt, re.S)
+    out = {}
+    for name, body in blocks:
+        d = {}
+        body_content = re.sub(r"!.*", "", body)  # remove comments
+        for line in body_content.split(","):
+            if "=" in line:
+                k, v = line.split("=", 1)
+                v = v.replace("T", "True").replace("F", "False")
+                try:
+                    v = ast.literal_eval(v.strip())
+                except Exception:
+                    v = v.strip().strip("'")
+                d[k.strip()] = v
+        out[name] = d
+    return out