diff --git a/model/atmosphere/dycore/tests/dycore/fixtures.py b/model/atmosphere/dycore/tests/dycore/fixtures.py index 6fc5afa77d..77c0f0c030 100644 --- a/model/atmosphere/dycore/tests/dycore/fixtures.py +++ b/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, ) diff --git a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py index f00a231235..d9bcfab3e4 100644 --- a/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py +++ b/model/atmosphere/dycore/tests/dycore/integration_tests/test_solve_nonhydro.py @@ -45,9 +45,9 @@ def test_validate_divdamp_fields_against_savepoint_values( grid_savepoint: sb.IconGridSavepoint, savepoint_nonhydro_init: sb.IconNonHydroInitSavepoint, icon_grid: base_grid.Grid, + solve_nonhydro_config: solve_nh.NonHydrostaticConfig, backend: gtx_typing.Backend, ) -> None: - config = solve_nh.NonHydrostaticConfig() second_order_divdamp_factor = 0.032 mean_cell_area = grid_savepoint.mean_cell_area() interpolated_fourth_order_divdamp_factor = data_alloc.zero_field( @@ -67,20 +67,20 @@ def test_validate_divdamp_fields_against_savepoint_values( ) smagorinsky.en_smag_fac_for_zero_nshift.with_backend(backend)( grid_savepoint.vct_a(), - config.fourth_order_divdamp_factor, - config.fourth_order_divdamp_factor2, - config.fourth_order_divdamp_factor3, - config.fourth_order_divdamp_factor4, - config.fourth_order_divdamp_z, - config.fourth_order_divdamp_z2, - config.fourth_order_divdamp_z3, - config.fourth_order_divdamp_z4, + solve_nonhydro_config.fourth_order_divdamp_factor, + solve_nonhydro_config.fourth_order_divdamp_factor2, + solve_nonhydro_config.fourth_order_divdamp_factor3, + solve_nonhydro_config.fourth_order_divdamp_factor4, + solve_nonhydro_config.fourth_order_divdamp_z, + solve_nonhydro_config.fourth_order_divdamp_z2, + solve_nonhydro_config.fourth_order_divdamp_z3, + solve_nonhydro_config.fourth_order_divdamp_z4, interpolated_fourth_order_divdamp_factor, offset_provider={"Koff": dims.KDim}, ) dycore_utils._calculate_fourth_order_divdamp_scaling_coeff.with_backend(backend)( interpolated_fourth_order_divdamp_factor=interpolated_fourth_order_divdamp_factor, - divdamp_order=config.divdamp_order, + divdamp_order=solve_nonhydro_config.divdamp_order, mean_cell_area=mean_cell_area, second_order_divdamp_factor=second_order_divdamp_factor, out=fourth_order_divdamp_scaling_coeff, @@ -90,7 +90,7 @@ def test_validate_divdamp_fields_against_savepoint_values( backend )( fourth_order_divdamp_scaling_coeff, - config.max_nudging_coefficient, + solve_nonhydro_config.max_nudging_coefficient, constants.DBL_EPS, out=reduced_fourth_order_divdamp_coeff_at_nest_boundary, offset_provider={}, @@ -172,24 +172,18 @@ def test_nonhydro_predictor_step( interpolation_savepoint, savepoint_nonhydro_exit, experiment, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, at_initial_timestep, caplog, backend, ): caplog.set_level(logging.WARN) - config = definitions.construct_nonhydrostatic_config(experiment) sp = savepoint_nonhydro_init sp_exit = savepoint_nonhydro_exit - nonhydro_params = solve_nh.NonHydrostaticParams(config) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + nonhydro_params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) dtime = sp.get_metadata("dtime").get("dtime") diagnostic_state_nh = utils.construct_diagnostics(sp, icon_grid, backend) @@ -202,7 +196,7 @@ def test_nonhydro_predictor_step( solve_nonhydro = solve_nh.SolveNonhydro( grid=icon_grid, - config=config, + config=solve_nonhydro_config, params=nonhydro_params, metric_state_nonhydro=metric_state_nonhydro, interpolation_state=interpolation_state, @@ -503,22 +497,16 @@ def test_nonhydro_corrector_step( interpolation_savepoint, savepoint_nonhydro_exit, experiment, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, caplog, backend, ): caplog.set_level(logging.WARN) - config = definitions.construct_nonhydrostatic_config(experiment) init_savepoint = savepoint_nonhydro_init - nonhydro_params = solve_nh.NonHydrostaticParams(config) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + nonhydro_params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) dtime = init_savepoint.get_metadata("dtime").get("dtime") lprep_adv = init_savepoint.get_metadata("prep_adv").get("prep_adv") prep_adv = dycore_states.PrepAdvection( @@ -552,7 +540,7 @@ def test_nonhydro_corrector_step( solve_nonhydro = solve_nh.SolveNonhydro( grid=icon_grid, - config=config, + config=solve_nonhydro_config, params=nonhydro_params, metric_state_nonhydro=metric_state_nonhydro, interpolation_state=interpolation_state, @@ -690,6 +678,8 @@ def test_run_solve_nonhydro_single_step( step_date_init, step_date_exit, experiment, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, icon_grid, savepoint_nonhydro_init, @@ -706,19 +696,11 @@ def test_run_solve_nonhydro_single_step( backend, ): caplog.set_level(logging.WARN) - config = definitions.construct_nonhydrostatic_config(experiment) sp = savepoint_nonhydro_init sp_step_exit = savepoint_nonhydro_step_final - nonhydro_params = solve_nh.NonHydrostaticParams(config) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + nonhydro_params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) dtime = sp.get_metadata("dtime").get("dtime") lprep_adv = sp.get_metadata("prep_adv").get("prep_adv") prep_adv = dycore_states.PrepAdvection( @@ -740,7 +722,7 @@ def test_run_solve_nonhydro_single_step( solve_nonhydro = solve_nh.SolveNonhydro( grid=icon_grid, - config=config, + config=solve_nonhydro_config, params=nonhydro_params, metric_state_nonhydro=metric_state_nonhydro, interpolation_state=interpolation_state, @@ -822,6 +804,8 @@ def test_run_solve_nonhydro_multi_step( at_initial_timestep, *, icon_grid, + vertical_grid_config, + solve_nonhydro_config, savepoint_nonhydro_init, lowest_layer_thickness, model_top_height, @@ -835,18 +819,10 @@ def test_run_solve_nonhydro_multi_step( ndyn_substeps, backend, ): - config = definitions.construct_nonhydrostatic_config(experiment) sp = savepoint_nonhydro_init sp_step_exit = savepoint_nonhydro_step_final - nonhydro_params = solve_nh.NonHydrostaticParams(config) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + nonhydro_params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) dtime = sp.get_metadata("dtime").get("dtime") lprep_adv = sp.get_metadata("prep_adv").get("prep_adv") prep_adv = dycore_states.PrepAdvection( @@ -873,7 +849,7 @@ def test_run_solve_nonhydro_multi_step( solve_nonhydro = solve_nh.SolveNonhydro( grid=icon_grid, - config=config, + config=solve_nonhydro_config, params=nonhydro_params, metric_state_nonhydro=metric_state_nonhydro, interpolation_state=interpolation_state, @@ -979,9 +955,8 @@ def test_run_solve_nonhydro_multi_step( @pytest.mark.datatest @pytest.mark.parametrize("experiment", [definitions.Experiments.MCH_CH_R04B09]) -def test_non_hydrostatic_params(savepoint_nonhydro_init): - config = solve_nh.NonHydrostaticConfig() - params = solve_nh.NonHydrostaticParams(config) +def test_non_hydrostatic_params(savepoint_nonhydro_init, solve_nonhydro_config): + params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) assert params.advection_implicit_weight_parameter == savepoint_nonhydro_init.wgt_nnew_vel() assert params.advection_explicit_weight_parameter == savepoint_nonhydro_init.wgt_nnow_vel() @@ -1013,6 +988,8 @@ def test_compute_perturbed_quantities_and_interpolation( step_date_init, step_date_exit, *, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, icon_grid, lowest_layer_thickness, @@ -1032,14 +1009,7 @@ def test_compute_perturbed_quantities_and_interpolation( sp_init = savepoint_nonhydro_init sp_ref = savepoint_compute_edge_diagnostics_for_dycore_and_update_vn_init sp_exit = savepoint_nonhydro_exit - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) current_rho = sp_init.rho_now() current_theta_v = sp_init.theta_v_now() @@ -1074,8 +1044,7 @@ def test_compute_perturbed_quantities_and_interpolation( icon_grid, dims.CellDim, dims.KDim, allocator=backend ) - config = definitions.construct_nonhydrostatic_config(experiment) - igradp_method = config.igradp_method + igradp_method = solve_nonhydro_config.igradp_method nflatlev = vertical_params.nflatlev nflat_gradp = grid_savepoint.nflat_gradp() @@ -1370,6 +1339,8 @@ def test_compute_rho_theta_pgrad_and_update_vn( step_date_init, step_date_exit, *, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, icon_grid, savepoint_nonhydro_init, @@ -1405,14 +1376,7 @@ def test_compute_rho_theta_pgrad_and_update_vn( end_edge_halo_level_2 = icon_grid.end_index(edge_domain(h_grid.Zone.HALO_LEVEL_2)) end_edge_local = icon_grid.end_index(edge_domain(h_grid.Zone.LOCAL)) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) current_vn = sp_stencil_init.vn() next_vn = sp_nh_init.vn_new() @@ -1438,15 +1402,14 @@ def test_compute_rho_theta_pgrad_and_update_vn( grf_tend_vn = sp_nh_init.grf_tend_vn() rho_at_edges_on_model_levels = sp_stencil_init.z_rho_e() theta_v_at_edges_on_model_levels = sp_stencil_init.z_theta_v_e() - config = definitions.construct_nonhydrostatic_config(experiment) primal_normal_cell_1 = grid_savepoint.primal_normal_cell_x() primal_normal_cell_2 = grid_savepoint.primal_normal_cell_y() dual_normal_cell_1 = grid_savepoint.dual_normal_cell_x() dual_normal_cell_2 = grid_savepoint.dual_normal_cell_y() - iau_wgt_dyn = config.iau_wgt_dyn - is_iau_active = config.is_iau_active - igradp_method = config.igradp_method + iau_wgt_dyn = solve_nonhydro_config.iau_wgt_dyn + is_iau_active = solve_nonhydro_config.is_iau_active + igradp_method = solve_nonhydro_config.igradp_method z_rho_e_ref = sp_stencil_exit.z_rho_e() z_theta_v_e_ref = sp_stencil_exit.z_theta_v_e() @@ -1586,6 +1549,7 @@ def test_apply_divergence_damping_and_update_vn( step_date_init, step_date_exit, *, + solve_nonhydro_config, ndyn_substeps, icon_grid, savepoint_nonhydro_init, @@ -1620,7 +1584,6 @@ def test_apply_divergence_damping_and_update_vn( current_vn = sp_stencil_init.vn() next_vn = savepoint_nonhydro_init.vn_new() horizontal_gradient_of_normal_wind_divergence = sp_nh_init.z_graddiv_vn() - config = definitions.construct_nonhydrostatic_config(experiment) mean_cell_area = grid_savepoint.mean_cell_area() # TODO: Use serialized data ('enh_divdamp_fac' in icon) instead of computing 'interpolated_fourth_order_divdamp_factor' @@ -1630,8 +1593,8 @@ def test_apply_divergence_damping_and_update_vn( allocator=backend, ) - iau_wgt_dyn = config.iau_wgt_dyn - divdamp_order = config.divdamp_order + iau_wgt_dyn = solve_nonhydro_config.iau_wgt_dyn + divdamp_order = solve_nonhydro_config.divdamp_order second_order_divdamp_scaling_coeff = sp_nh_init.divdamp_fac_o2() * mean_cell_area second_order_divdamp_factor = savepoint_nonhydro_init.divdamp_fac_o2() apply_2nd_order_divergence_damping = ( @@ -1642,23 +1605,25 @@ def test_apply_divergence_damping_and_update_vn( divdamp_order == dycore_states.DivergenceDampingOrder.FOURTH_ORDER or ( divdamp_order == dycore_states.DivergenceDampingOrder.COMBINED - and second_order_divdamp_factor <= (4.0 * config.fourth_order_divdamp_factor) + and second_order_divdamp_factor + <= (4.0 * solve_nonhydro_config.fourth_order_divdamp_factor) ) ) - is_iau_active = config.is_iau_active + # TODO (Chia Rui): a bug to be fixed in https://github.com/C2SM/icon4py/pull/972, is_iau_active is a runtime variable and it is not serialized + is_iau_active = False vn_ref = sp_nh_exit.vn_new() smagorinsky.en_smag_fac_for_zero_nshift.with_backend(backend)( grid_savepoint.vct_a(), - config.fourth_order_divdamp_factor, - config.fourth_order_divdamp_factor2, - config.fourth_order_divdamp_factor3, - config.fourth_order_divdamp_factor4, - config.fourth_order_divdamp_z, - config.fourth_order_divdamp_z2, - config.fourth_order_divdamp_z3, - config.fourth_order_divdamp_z4, + solve_nonhydro_config.fourth_order_divdamp_factor, + solve_nonhydro_config.fourth_order_divdamp_factor2, + solve_nonhydro_config.fourth_order_divdamp_factor3, + solve_nonhydro_config.fourth_order_divdamp_factor4, + solve_nonhydro_config.fourth_order_divdamp_z, + solve_nonhydro_config.fourth_order_divdamp_z2, + solve_nonhydro_config.fourth_order_divdamp_z3, + solve_nonhydro_config.fourth_order_divdamp_z4, interpolated_fourth_order_divdamp_factor, offset_provider={"Koff": dims.KDim}, ) @@ -1694,7 +1659,7 @@ def test_apply_divergence_damping_and_update_vn( divdamp_order=divdamp_order, mean_cell_area=mean_cell_area, second_order_divdamp_factor=second_order_divdamp_factor, - max_nudging_coefficient=config.max_nudging_coefficient, + max_nudging_coefficient=solve_nonhydro_config.max_nudging_coefficient, dbl_eps=constants.DBL_EPS, horizontal_start=start_edge_nudging_level_2, horizontal_end=end_edge_local, @@ -1733,14 +1698,13 @@ def test_apply_divergence_damping_and_update_vn( ], ) def test_compute_horizontal_velocity_quantities_and_fluxes( - istep_init, - istep_exit, - substep_init, - substep_exit, + experiment, step_date_init, step_date_exit, - experiment, + *, icon_grid, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, grid_savepoint, lowest_layer_thickness, @@ -1756,14 +1720,7 @@ def test_compute_horizontal_velocity_quantities_and_fluxes( backend, ): edge_domain = h_grid.domain(dims.EdgeDim) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) ddqz_z_full_e = metrics_savepoint.ddqz_z_full_e() ddxn_z_full = metrics_savepoint.ddxn_z_full() @@ -2036,6 +1993,8 @@ def test_vertically_implicit_solver_at_predictor_step( step_date_init, step_date_exit, *, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, icon_grid, savepoint_nonhydro_init, @@ -2055,17 +2014,9 @@ def test_vertically_implicit_solver_at_predictor_step( ): sp_nh_exit = savepoint_nonhydro_exit sp_stencil_init = savepoint_vertically_implicit_dycore_solver_init - config = definitions.construct_nonhydrostatic_config(experiment) xp = data_alloc.import_array_ns(backend) - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) at_first_substep = substep_init == 1 contravariant_correction_at_edges_on_model_levels = sp_nh_exit.z_w_concorr_me() @@ -2092,9 +2043,10 @@ def test_vertically_implicit_solver_at_predictor_step( dwdz_at_cells_on_model_levels = sp_stencil_init.z_dwdz_dd() exner_dynamical_increment = sp_stencil_init.exner_dyn_incr() - iau_wgt_dyn = config.iau_wgt_dyn - is_iau_active = config.is_iau_active - divdamp_type = config.divdamp_type + # TODO (Chia Rui): a bug to be fixed in https://github.com/C2SM/icon4py/pull/972, is_iau_active is a runtime variable and it is not serialized + iau_wgt_dyn = 0.0 + is_iau_active = False + divdamp_type = solve_nonhydro_config.divdamp_type w_concorr_c_ref = sp_nh_exit.w_concorr_c() w_ref = sp_nh_exit.w_new() @@ -2157,7 +2109,7 @@ def test_vertically_implicit_solver_at_predictor_step( iau_wgt_dyn=iau_wgt_dyn, dtime=savepoint_nonhydro_init.get_metadata("dtime").get("dtime"), is_iau_active=is_iau_active, - rayleigh_type=config.rayleigh_type, + rayleigh_type=solve_nonhydro_config.rayleigh_type, divdamp_type=divdamp_type, at_first_substep=at_first_substep, end_index_of_damping_layer=grid_savepoint.nrdmax(), @@ -2197,7 +2149,7 @@ def test_vertically_implicit_solver_at_predictor_step( # manually set the reference equal to zero when k < starting_vertical_index_for_3d_divdamp. starting_vertical_index_for_3d_divdamp = ( xp.min(xp.where(metrics_savepoint.scaling_factor_for_3d_divdamp().ndarray > 0.0))[0] - if config.divdamp_type == 32 + if solve_nonhydro_config.divdamp_type == 32 else 0 ) z_dwdz_dd_ref_with_zero_in_2d_divdamp_layers = z_dwdz_dd_ref.asnumpy() @@ -2241,6 +2193,8 @@ def test_vertically_implicit_solver_at_corrector_step( step_date_init, step_date_exit, *, + vertical_grid_config, + solve_nonhydro_config, ndyn_substeps, icon_grid, savepoint_nonhydro_init, @@ -2257,20 +2211,12 @@ def test_vertically_implicit_solver_at_corrector_step( ): sp_nh_exit = savepoint_nonhydro_exit sp_stencil_init = savepoint_vertically_implicit_dycore_solver_init - vertical_config = v_grid.VerticalGridConfig( - icon_grid.num_levels, - lowest_layer_thickness=lowest_layer_thickness, - model_top_height=model_top_height, - stretch_factor=stretch_factor, - rayleigh_damping_height=damping_height, - ) - vertical_params = utils.create_vertical_params(vertical_config, grid_savepoint) + vertical_params = utils.create_vertical_params(vertical_grid_config, grid_savepoint) at_first_substep = substep_init == 0 at_last_substep = substep_exit == 0 - config = definitions.construct_nonhydrostatic_config(experiment) - nonhydro_params = solve_nh.NonHydrostaticParams(config) + nonhydro_params = solve_nh.NonHydrostaticParams(solve_nonhydro_config) mass_flux_at_edges_on_model_levels = sp_stencil_init.mass_fl_e() theta_v_flux_at_edges_on_model_levels = sp_stencil_init.z_theta_v_fl_e() @@ -2301,8 +2247,9 @@ def test_vertically_implicit_solver_at_corrector_step( r_nsubsteps = 1.0 / ndyn_substeps kstart_moist = vertical_params.kstart_moist - iau_wgt_dyn = config.iau_wgt_dyn - is_iau_active = config.is_iau_active + # TODO (Chia Rui): a bug to be fixed in https://github.com/C2SM/icon4py/pull/972, is_iau_active is a runtime variable and it is not serialized + iau_wgt_dyn = 0.0 + is_iau_active = False w_ref = sp_nh_exit.w_new() rho_ref = sp_nh_exit.rho_new() @@ -2364,7 +2311,7 @@ def test_vertically_implicit_solver_at_corrector_step( iau_wgt_dyn=iau_wgt_dyn, dtime=savepoint_nonhydro_init.get_metadata("dtime").get("dtime"), is_iau_active=is_iau_active, - rayleigh_type=config.rayleigh_type, + rayleigh_type=solve_nonhydro_config.rayleigh_type, at_first_substep=at_first_substep, at_last_substep=at_last_substep, end_index_of_damping_layer=grid_savepoint.nrdmax(), diff --git a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/microphysics_options.py b/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/microphysics_options.py index b41b2777cd..48b5bce8a3 100644 --- a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/microphysics_options.py +++ b/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, +] diff --git a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/single_moment_six_class_gscp_graupel.py b/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/single_moment_six_class_gscp_graupel.py index a62413314a..20f9863de7 100644 --- a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/single_moment_six_class_gscp_graupel.py +++ b/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/single_moment_six_class_gscp_graupel.py @@ -61,11 +61,11 @@ 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. @@ -73,7 +73,7 @@ class SingleMomentSixClassIconGraupelConfig: #: 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. diff --git a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/stencils/microphyiscal_processes.py b/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/stencils/microphyiscal_processes.py index 7948cfc940..c4aa43438e 100644 --- a/model/atmosphere/subgrid_scale_physics/microphysics/src/icon4py/model/atmosphere/subgrid_scale_physics/microphysics/stencils/microphyiscal_processes.py +++ b/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 diff --git a/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_saturation_adjustment.py b/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_saturation_adjustment.py index 12af8065ae..108ee0e665 100644 --- a/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_saturation_adjustment.py +++ b/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, diff --git a/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_single_moment_six_class_gscp_graupel.py b/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_single_moment_six_class_gscp_graupel.py index 6c571a0e54..b2d03922e4 100644 --- a/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_single_moment_six_class_gscp_graupel.py +++ b/model/atmosphere/subgrid_scale_physics/microphysics/tests/microphysics/integration_tests/test_single_moment_six_class_gscp_graupel.py @@ -38,9 +38,9 @@ @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( @@ -48,20 +48,17 @@ ) 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, diff --git a/model/testing/pyproject.toml b/model/testing/pyproject.toml index 9d7b94d752..cdf5b59759 100644 --- a/model/testing/pyproject.toml +++ b/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" diff --git a/model/testing/src/icon4py/model/testing/datatest_utils.py b/model/testing/src/icon4py/model/testing/datatest_utils.py index b8741cbaa1..8b01e5d86b 100644 --- a/model/testing/src/icon4py/model/testing/datatest_utils.py +++ b/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 diff --git a/model/testing/src/icon4py/model/testing/definitions.py b/model/testing/src/icon4py/model/testing/definitions.py index 0668ffca5a..66b01bc823 100644 --- a/model/testing/src/icon4py/model/testing/definitions.py +++ b/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 _: diff --git a/model/testing/src/icon4py/model/testing/fixtures/datatest.py b/model/testing/src/icon4py/model/testing/fixtures/datatest.py index 5d59060a33..91be067439 100644 --- a/model/testing/src/icon4py/model/testing/fixtures/datatest.py +++ b/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 diff --git a/uv.lock b/uv.lock index 4b38bfa69d..91e669dc90 100644 --- a/uv.lock +++ b/uv.lock @@ -1152,6 +1152,16 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/b5/fd/afcd0496feca3276f509df3dbd5dae726fcc756f1a08d9e25abe1733f962/executing-2.1.0-py2.py3-none-any.whl", hash = "sha256:8d63781349375b5ebccc3142f4b30350c0cd9c79f921cde38be2be4637e98eaf", size = 25805 }, ] +[[package]] +name = "f90nml" +version = "1.5" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/54/39/5bef20ee2868924029b2d34523839798cc8e48183ca57d9704ae6ef92da6/f90nml-1.5.tar.gz", hash = "sha256:d5585f0051234c742f7515c4b1f8484f8f628a71c4e31f2602357089dcc05fc1", size = 70959, upload_time = "2025-10-07T15:07:14.365Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/a9/ee/39fc4b67a7ae5ba91d6da18a1b66cbdf6b2f4f88e3a26ca36823f3482dd3/f90nml-1.5-py2.py3-none-any.whl", hash = "sha256:38bbbdfe81eb7ac62e924021c4df6bb5cf4846e2f8ad79bfecb0fd654250180a", size = 51971, upload_time = "2025-10-07T15:07:13.332Z" }, + { url = "https://files.pythonhosted.org/packages/21/f2/4454eefc15cc326b46530d230c58cc0bb91a1e9797f2842b2a1720cbb233/f90nml-1.5.0-py2.py3-none-any.whl", hash = "sha256:bdf616dbe7e83619feb86d54358fb8d97038133bfd8f9ba9a01eeca5dc4691a7", size = 51994, upload_time = "2025-10-07T15:25:09.064Z" }, +] + [[package]] name = "factory-boy" version = "3.3.3" @@ -1962,6 +1972,7 @@ name = "icon4py-testing" version = "0.0.6" source = { editable = "model/testing" } dependencies = [ + { name = "f90nml" }, { name = "filelock" }, { name = "gt4py" }, { name = "icon4py-common", extra = ["io"] }, @@ -1975,6 +1986,7 @@ dependencies = [ [package.metadata] requires-dist = [ + { name = "f90nml", specifier = ">=1.5" }, { name = "filelock", specifier = ">=3.18.0" }, { name = "gt4py", specifier = "==1.1.4" }, { name = "icon4py-common", extras = ["io"], editable = "model/common" },