Skip to content

Commit

Permalink
remove redundant code
Browse files Browse the repository at this point in the history
  • Loading branch information
OngChia committed Sep 4, 2024
1 parent 8a3f94a commit 73655e8
Show file tree
Hide file tree
Showing 6 changed files with 93 additions and 103 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,11 @@
where,
)

from icon4py.model.atmosphere.subgrid_scale_physics.microphysics import saturation_adjustment as satad
from icon4py.model.atmosphere.subgrid_scale_physics.microphysics import (
saturation_adjustment as satad,
)
from icon4py.model.common.dimension import CellDim, KDim
from icon4py.model.common.grid import icon as icon_grid, vertical as v_grid, horizontal as h_grid
from icon4py.model.common.grid import horizontal as h_grid, icon as icon_grid, vertical as v_grid
from icon4py.model.common.settings import backend
from icon4py.model.common.states import (
diagnostic_state as diagnostics,
Expand Down Expand Up @@ -257,13 +259,13 @@ class SingleMomentSixClassIconGraupelParams(FrozenNamespace):
latent_heat_for_vaporisation: wpfloat = 2.5008e6
#: Latent heat of sublimation for water [J/kg]. Originally expressed as als in ICON.
latent_heat_for_sublimation: wpfloat = 2.8345e6
#: Melting temperature of ice/snow [K]
tmelt: wpfloat = 273.15
#: Melting temperature of ice/snow [K]. Originally expressed as tmelt in ICON.
melting_temperature: wpfloat = 273.15
#: Triple point of water at 611hPa [K]
t3: wpfloat = 273.16

#: ice crystal number concentration at threshold temperature for mixed-phase cloud
nimix: wpfloat = 5.0 * np.exp(0.304 * (tmelt - threshold_freeze_temperature_mixedphase))
nimix: wpfloat = 5.0 * np.exp(0.304 * (melting_temperature - threshold_freeze_temperature_mixedphase))

#: Gas constant of dry air [J/K/kg]
rd: wpfloat = 287.04
Expand All @@ -290,7 +292,7 @@ class SingleMomentSixClassIconGraupelParams(FrozenNamespace):
ccshi1: wpfloat = (
latent_heat_for_sublimation * latent_heat_for_sublimation / (dry_air_latent_heat * rv)
)
ccdvtp: wpfloat = 2.22e-5 * tmelt ** (-1.94) * 101325.0
ccdvtp: wpfloat = 2.22e-5 * melting_temperature ** (-1.94) * 101325.0
ccidep: wpfloat = 4.0 * ice_m0 ** (-1.0 / 3.0)
ccswxp_ln1o2: wpfloat = np.exp(ccswxp * np.log(0.5))

Expand All @@ -306,15 +308,15 @@ class SingleMomentSixClassIconGraupelParams(FrozenNamespace):
#: [K*kg/J]"""
rcvd: wpfloat = 1.0 / cvd

pvsw0: wpfloat = tetens_p0 * np.exp(tetens_aw * (tmelt - tmelt) / (tmelt - tetens_bw))
pvsw0: wpfloat = tetens_p0 * np.exp(tetens_aw * (melting_temperature - melting_temperature) / (melting_temperature - tetens_bw))


icon_graupel_params: Final = SingleMomentSixClassIconGraupelParams()


@gtx.field_operator
def _compute_cooper_inp_concentration(temperature: wpfloat) -> wpfloat:
cnin = 5.0 * exp(0.304 * (icon_graupel_params.tmelt - temperature))
cnin = 5.0 * exp(0.304 * (icon_graupel_params.melting_temperature - temperature))
cnin = minimum(cnin, icon_graupel_params.nimax_Thom)
return cnin

Expand Down Expand Up @@ -357,17 +359,19 @@ def _compute_snow_interception_and_collision_parameters(
if snow_intercept_option == 1:
# Calculate n0s using the temperature-dependent
# formula of Field et al. (2005)
local_tc = temperature - icon_graupel_params.tmelt
local_tc = temperature - icon_graupel_params.melting_temperature
local_tc = minimum(local_tc, wpfloat("0.0"))
local_tc = maximum(local_tc, wpfloat("-40.0"))
n0s = icon_graupel_params.snow_intercept_parameter_n0s1 * exp(icon_graupel_params.snow_intercept_parameter_n0s2 * local_tc)
n0s = icon_graupel_params.snow_intercept_parameter_n0s1 * exp(
icon_graupel_params.snow_intercept_parameter_n0s2 * local_tc
)
n0s = minimum(n0s, wpfloat("1.0e9"))
n0s = maximum(n0s, wpfloat("1.0e6"))

elif snow_intercept_option == 2:
# Calculate n0s using the temperature-dependent moment
# relations of Field et al. (2005)
local_tc = temperature - icon_graupel_params.tmelt
local_tc = temperature - icon_graupel_params.melting_temperature
local_tc = minimum(local_tc, wpfloat("0.0"))
local_tc = maximum(local_tc, wpfloat("-40.0"))

Expand Down Expand Up @@ -662,7 +666,7 @@ def _riming_in_clouds(

srimg_c2g = icon_graupel_params.crim_g * qc * celnrimexp_g

if temperature >= icon_graupel_params.tmelt:
if temperature >= icon_graupel_params.melting_temperature:
sshed_c2r = srims_c2s + srimg_c2g
srims_c2s = wpfloat("0.0")
srimg_c2g = wpfloat("0.0")
Expand Down Expand Up @@ -829,11 +833,11 @@ def _collision_and_ice_deposition_in_cold_ice_clouds(
Returns:
cloud-to-rain autoconversionn rate, rain-cloud accretion rate
"""
if (temperature <= icon_graupel_params.tmelt) & llqi:
if (temperature <= icon_graupel_params.melting_temperature) & llqi:
# Change in sticking efficiency needed in case of cloud ice sedimentation
# (based on Guenther Zaengls work)
local_eff = minimum(
exp(wpfloat("0.09") * (temperature - icon_graupel_params.tmelt)), wpfloat("1.0")
exp(wpfloat("0.09") * (temperature - icon_graupel_params.melting_temperature)), wpfloat("1.0")
)
local_eff = maximum(local_eff, ice_stickeff_min)
local_eff = maximum(
Expand Down Expand Up @@ -967,7 +971,7 @@ def _snow_and_graupel_depositional_growth_in_cold_ice_clouds(
depositional growth rate of snow, depositional growth rate of graupel
"""
if llqi | llqs | llqg:
if temperature <= icon_graupel_params.tmelt:
if temperature <= icon_graupel_params.melting_temperature:
local_qvsidiff = qv - qvsi
local_svmax = local_qvsidiff / dt

Expand Down Expand Up @@ -1054,24 +1058,24 @@ def _melting(
depositional growth rate of snow, depositional growth rate of graupel
"""
if llqi | llqs | llqg:
if temperature > icon_graupel_params.tmelt:
if temperature > icon_graupel_params.melting_temperature:
# cloud ice melts instantaneously
simlt_i2c = rhoqi_intermediate / rho / dt

local_qvsw0 = icon_graupel_params.pvsw0 / (
rho * icon_graupel_params.rv * icon_graupel_params.tmelt
rho * icon_graupel_params.rv * icon_graupel_params.melting_temperature
)
local_qvsw0diff = qv - local_qvsw0

# ** GZ: several numerical fits in this section should be replaced with physically more meaningful formulations **
if (
temperature
> icon_graupel_params.tmelt - icon_graupel_params.tcrit * local_qvsw0diff
> icon_graupel_params.melting_temperature - icon_graupel_params.tcrit * local_qvsw0diff
):
# calculate melting rate
local_x1 = (
temperature
- icon_graupel_params.tmelt
- icon_graupel_params.melting_temperature
+ icon_graupel_params.asmel * local_qvsw0diff
)
ssmlt_s2r = (
Expand Down Expand Up @@ -1200,7 +1204,7 @@ def _evaporation_and_freezing_in_subsaturated_air(
local_lnqr = log(rhoqr)
local_x1 = wpfloat("1.0") + bev * exp(bevxp * local_lnqr)
# Limit evaporation rate in order to avoid overshoots towards supersaturation, the pre-factor approximates (esat(T_wb)-e)/(esat(T)-e) at temperatures between 0 degC and 30 degC
local_temp_c = temperature - icon_graupel_params.tmelt
local_temp_c = temperature - icon_graupel_params.melting_temperature
local_maxevap = (
(
wpfloat("0.61")
Expand Down Expand Up @@ -1241,7 +1245,7 @@ def _evaporation_and_freezing_in_subsaturated_air(
def sat_pres_water(temperature: wpfloat) -> wpfloat:
return icon_graupel_params.tetens_p0 * exp(
icon_graupel_params.tetens_aw
* (temperature - icon_graupel_params.tmelt)
* (temperature - icon_graupel_params.melting_temperature)
/ (temperature - icon_graupel_params.tetens_bw)
)

Expand All @@ -1250,7 +1254,7 @@ def sat_pres_water(temperature: wpfloat) -> wpfloat:
def sat_pres_ice(temperature: wpfloat) -> wpfloat:
return icon_graupel_params.tetens_p0 * exp(
icon_graupel_params.tetens_ai
* (temperature - icon_graupel_params.tmelt)
* (temperature - icon_graupel_params.melting_temperature)
/ (temperature - icon_graupel_params.tetens_bi)
)

Expand All @@ -1275,7 +1279,12 @@ def __init__(
self.metric_state = metric_state
self.vertical_params = vertical_params
self.saturation_adjustment = satad.SaturationAdjustment(
config=saturation_adjust_config, grid=grid, vertical_params=vertical_params,metric_state=satad.MetricStateSaturationAdjustment(ddqz_z_full=metric_state.ddqz_z_full),
config=saturation_adjust_config,
grid=grid,
vertical_params=vertical_params,
metric_state=satad.MetricStateSaturationAdjustment(
ddqz_z_full=metric_state.ddqz_z_full
),
)

self._initialize_local_fields()
Expand Down Expand Up @@ -2006,15 +2015,15 @@ def _icon_graupel_scan(
if use_constant_water_heat_capacity
else icon_graupel_params.latent_heat_for_vaporisation
+ (icon_graupel_params.cp_v - icon_graupel_params.spec_heat_cap_water)
* (temperature - icon_graupel_params.tmelt)
* (temperature - icon_graupel_params.melting_temperature)
- icon_graupel_params.rv * temperature
)
lhs = (
icon_graupel_params.latent_heat_for_sublimation
if use_constant_water_heat_capacity
else icon_graupel_params.latent_heat_for_sublimation
+ (icon_graupel_params.cp_v - icon_graupel_params.spec_heat_cap_ice)
* (temperature - icon_graupel_params.tmelt)
* (temperature - icon_graupel_params.melting_temperature)
- icon_graupel_params.rv * temperature
)

Expand Down Expand Up @@ -2410,7 +2419,7 @@ def _icon_graupel_scan(
scosg_s2g = minimum(scosg_s2g, srims_c2s + cssmax)

if llqi | llqs | llqg:
if temperature <= icon_graupel_params.tmelt: # cold case
if temperature <= icon_graupel_params.melting_temperature: # cold case
qvsidiff = qv - qvsi
csimax = rhoqi_intermediate / rho / dt

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -122,40 +122,51 @@ def test_graupel(
assert helpers.dallclose(graupel_microphysics.ccs[0], init_savepoint.ccsrim(), atol=1.0e-8)
assert helpers.dallclose(graupel_microphysics.ccs[1], init_savepoint.ccsagg(), atol=1.0e-8)
assert helpers.dallclose(graupel.icon_graupel_params.ccsaxp, init_savepoint.ccsaxp())
assert helpers.dallclose(graupel.icon_graupel_params.ccsdep, init_savepoint.ccsdep(), atol=1.0e-7)
assert helpers.dallclose(
graupel.icon_graupel_params.ccsdep, init_savepoint.ccsdep(), atol=1.0e-7
)
assert helpers.dallclose(graupel.icon_graupel_params.ccsdxp, init_savepoint.ccsdxp())
assert helpers.dallclose(graupel_microphysics.ccs[2], init_savepoint.ccsvel(), atol=1.0e-8)
assert helpers.dallclose(graupel.icon_graupel_params.ccsvxp, init_savepoint.ccsvxp())
assert helpers.dallclose(graupel.icon_graupel_params.ccslam, init_savepoint.ccslam(), atol=1.0e-10)
assert helpers.dallclose(
graupel.icon_graupel_params.ccslam, init_savepoint.ccslam(), atol=1.0e-10
)
assert helpers.dallclose(graupel.icon_graupel_params.ccslxp, init_savepoint.ccslxp())
assert helpers.dallclose(graupel.icon_graupel_params.ccshi1, init_savepoint.ccshi1())
assert helpers.dallclose(graupel.icon_graupel_params.ccdvtp, init_savepoint.ccdvtp())
assert helpers.dallclose(graupel.icon_graupel_params.ccidep, init_savepoint.ccidep())
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[0], init_savepoint.vzxp())
assert helpers.dallclose(
graupel_microphysics.rain_vel_coef[1], init_savepoint.vz0r(), atol=1.0e-10
)
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[2], init_savepoint.cevxp())
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[3], init_savepoint.cev(), atol=1.0e-10)
assert helpers.dallclose(
graupel_microphysics.rain_vel_coef[3], init_savepoint.cev(), atol=1.0e-10
)
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[4], init_savepoint.bevxp())
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[5], init_savepoint.bev())
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[0], init_savepoint.vzxp())
assert helpers.dallclose(graupel_microphysics.rain_vel_coef[1], init_savepoint.vz0r(), atol=1.0e-10)

# TODO (Chia Rui): remove this slicing process, which finds the column with maximum tendency, when either the scan operator can be run on the gtfn backend or running on embedded backend is faster
max_index = np.unravel_index(
np.abs(tracer_state.qv.ndarray - exit_savepoint.qv().ndarray).argmax(),
exit_savepoint.qv().ndarray.shape,
)
cell_lower_limit = max_index[0] - 300 # - 277
cell_upper_limit = max_index[0] + 300 # - 276
cell_lower_limit = max_index[0] - 300
cell_upper_limit = max_index[0] + 300
cell_size = cell_upper_limit - cell_lower_limit
vertical_size = icon_grid.num_levels

slice_t = gtx.as_field(
(dims.CellDim, dims.KDim), diagnostic_state.temperature.ndarray[cell_lower_limit:cell_upper_limit, :]
(dims.CellDim, dims.KDim),
diagnostic_state.temperature.ndarray[cell_lower_limit:cell_upper_limit, :],
)
slice_pres = gtx.as_field(
(dims.CellDim, dims.KDim), diagnostic_state.pressure.ndarray[cell_lower_limit:cell_upper_limit, :]
(dims.CellDim, dims.KDim),
diagnostic_state.pressure.ndarray[cell_lower_limit:cell_upper_limit, :],
)
slice_rho = gtx.as_field(
(dims.CellDim, dims.KDim), prognostic_state.rho.ndarray[cell_lower_limit:cell_upper_limit, :]
(dims.CellDim, dims.KDim),
prognostic_state.rho.ndarray[cell_lower_limit:cell_upper_limit, :],
)
slice_qv = gtx.as_field(
(dims.CellDim, dims.KDim), tracer_state.qv.ndarray[cell_lower_limit:cell_upper_limit, :]
Expand Down Expand Up @@ -183,13 +194,27 @@ def test_graupel(
graupel_microphysics.metric_state.ddqz_z_full.ndarray[cell_lower_limit:cell_upper_limit, :],
)

slice_t_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qv_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qc_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qi_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qr_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qs_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_qg_t = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_t_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qv_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qc_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qi_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qr_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qs_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_qg_t = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)

slice_rhoqrv_old_kup = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
Expand All @@ -203,10 +228,18 @@ def test_graupel(
slice_rhoqiv_old_kup = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_vnew_r = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_vnew_s = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_vnew_g = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_vnew_i = gtx.as_field((dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float))
slice_vnew_r = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_vnew_s = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_vnew_g = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_vnew_i = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
slice_rain_precipitation_flux = gtx.as_field(
(dims.CellDim, dims.KDim), np.zeros((cell_size, vertical_size), dtype=float)
)
Expand Down
3 changes: 0 additions & 3 deletions model/common/src/icon4py/model/common/constants.py
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,6 @@
GAS_CONSTANT_WATER_VAPOR: Final[wpfloat] = 461.51
RV: Final[wpfloat] = GAS_CONSTANT_WATER_VAPOR

#: Melting temperature of ice/snow [K]
TMELT: Final[wpfloat] = 273.15

#: RV/RD - 1, tvmpc1 in ICON.
RV_O_RD_MINUS_1: Final[wpfloat] = GAS_CONSTANT_WATER_VAPOR / GAS_CONSTANT_DRY_AIR - 1.0
TVMPC1: Final[wpfloat] = RV_O_RD_MINUS_1
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -85,11 +85,9 @@ def get_global_grid_params(experiment: str) -> tuple[int, int]:
root, level = map(int, re.search("[Rr](\d+)[Bb](\d+)", experiment).groups())
return root, level
except AttributeError as err:
#raise ValueError(
# f"Could not parse grid_root and grid_level from experiment: {experiment} no 'rXbY'pattern."
#) from err
# TODO (Chia Rui): remove this temporary fix after merging from Jacopo's branch
return 0, 2
raise ValueError(
f"Could not parse grid_root and grid_level from experiment: {experiment} no 'rXbY'pattern."
) from err


def get_grid_id_for_experiment(experiment) -> uuid.UUID:
Expand Down
Loading

0 comments on commit 73655e8

Please sign in to comment.