Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -113,12 +113,11 @@ def test_diffusion_benchmark(
stretch_factor=1.0,
rayleigh_damping_height=1.0,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, allocator=allocator)
vct_a = v_grid.get_vct_a(vertical_config, allocator=allocator)

vertical_grid = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)

interpolation_state = diffusion_states.DiffusionInterpolationState(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -194,11 +194,10 @@ def test_diffusion_init(
stretch_factor=stretch_factor,
rayleigh_damping_height=damping_height,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, backend)
vct_a = v_grid.get_vct_a(vertical_config, backend)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)

meta = savepoint_diffusion_init.get_metadata("linit", "date")
Expand Down Expand Up @@ -327,11 +326,10 @@ def test_verify_diffusion_init_against_savepoint(
stretch_factor=stretch_factor,
rayleigh_damping_height=damping_height,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, backend)
vct_a = v_grid.get_vct_a(vertical_config, backend)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)

diffusion_granule = diffusion.Diffusion(
Expand Down Expand Up @@ -411,11 +409,10 @@ def test_run_diffusion_single_step(
stretch_factor=stretch_factor,
rayleigh_damping_height=damping_height,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, backend)
vct_a = v_grid.get_vct_a(vertical_config, backend)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)

config = definitions.construct_diffusion_config(experiment, ndyn_substeps)
Expand Down Expand Up @@ -490,9 +487,7 @@ def test_run_diffusion_multiple_steps(
rayleigh_damping_height=damping_height,
)

vertical_params = v_grid.VerticalGrid(
config=vertical_config, vct_a=grid_savepoint.vct_a(), vct_b=grid_savepoint.vct_b()
)
vertical_params = v_grid.VerticalGrid(config=vertical_config, vct_a=grid_savepoint.vct_a())

config = definitions.construct_diffusion_config(experiment, ndyn_substeps)
additional_parameters = diffusion.DiffusionParams(config)
Expand Down Expand Up @@ -606,12 +601,8 @@ def test_run_diffusion_initial_step(
stretch_factor=stretch_factor,
rayleigh_damping_height=damping_height,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, backend)
vertical_grid = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)
vct_a = v_grid.get_vct_a(vertical_config, backend)
vertical_grid = v_grid.VerticalGrid(config=vertical_config, vct_a=vct_a)
diagnostic_state = diffusion_states.DiffusionDiagnosticState(
hdef_ic=savepoint_diffusion_init.hdef_ic(),
div_ic=savepoint_diffusion_init.div_ic(),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -99,11 +99,7 @@ def test_parallel_diffusion(
grid=icon_grid,
config=config,
params=diffusion_params,
vertical_grid=v_grid.VerticalGrid(
vertical_config,
grid_savepoint.vct_a(),
grid_savepoint.vct_b(),
),
vertical_grid=v_grid.VerticalGrid(vertical_config, grid_savepoint.vct_a()),
metric_state=metric_state,
interpolation_state=interpolation_state,
edge_params=edge_geometry,
Expand Down Expand Up @@ -231,11 +227,7 @@ def test_parallel_diffusion_multiple_steps(
grid=icon_grid,
config=config,
params=diffusion_params,
vertical_grid=v_grid.VerticalGrid(
vertical_config,
grid_savepoint.vct_a(),
grid_savepoint.vct_b(),
),
vertical_grid=v_grid.VerticalGrid(vertical_config, grid_savepoint.vct_a()),
metric_state=metric_state,
interpolation_state=interpolation_state,
edge_params=edge_geometry,
Expand Down Expand Up @@ -283,11 +275,7 @@ def test_parallel_diffusion_multiple_steps(
grid=icon_grid,
config=config,
params=diffusion_params,
vertical_grid=v_grid.VerticalGrid(
vertical_config,
grid_savepoint.vct_a(),
grid_savepoint.vct_b(),
),
vertical_grid=v_grid.VerticalGrid(vertical_config, grid_savepoint.vct_a()),
metric_state=metric_state,
interpolation_state=interpolation_state,
edge_params=edge_geometry,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,9 @@ def solve_nonhydro(
stretch_factor=1.0,
rayleigh_damping_height=1.0,
)
vct_a, vct_b = v_grid.get_vct_a_and_vct_b(vertical_config, allocator=allocator)
vct_a = v_grid.get_vct_a(vertical_config, allocator=allocator)

vertical_grid = v_grid.VerticalGrid(
config=vertical_config,
vct_a=vct_a,
vct_b=vct_b,
)
vertical_grid = v_grid.VerticalGrid(config=vertical_config, vct_a=vct_a)

cell_geometry = grid_states.CellParams(
cell_center_lat=geometry_field_source.get(geometry_meta.CELL_LAT),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,7 @@ def _compare_cfl(
def create_vertical_params(
vertical_config: v_grid.VerticalGridConfig, grid_savepoint: serialbox.IconGridSavepoint
) -> v_grid.VerticalGrid:
return v_grid.VerticalGrid(
config=vertical_config, vct_a=grid_savepoint.vct_a(), vct_b=grid_savepoint.vct_b()
)
return v_grid.VerticalGrid(config=vertical_config, vct_a=grid_savepoint.vct_a())


@pytest.mark.embedded_static_args
Expand Down
6 changes: 1 addition & 5 deletions model/atmosphere/dycore/tests/dycore/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -86,11 +86,7 @@ def create_vertical_params(
vertical_config: v_grid.VerticalGridConfig,
sp: sb.IconGridSavepoint,
) -> v_grid.VerticalGrid:
return v_grid.VerticalGrid(
config=vertical_config,
vct_a=sp.vct_a(),
vct_b=sp.vct_b(),
)
return v_grid.VerticalGrid(config=vertical_config, vct_a=sp.vct_a())


def construct_diagnostics(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,7 @@ def test_saturation_adjustement(
icon_grid.num_levels,
model_top_height=model_top_height,
)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=grid_savepoint.vct_a(),
vct_b=grid_savepoint.vct_b(),
)
vertical_params = v_grid.VerticalGrid(config=vertical_config, vct_a=grid_savepoint.vct_a())
temperature_tendency = data_alloc.zero_field(
icon_grid, dims.CellDim, dims.KDim, allocator=backend
)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -64,11 +64,7 @@ def test_graupel(
lowest_layer_thickness=lowest_layer_thickness,
model_top_height=model_top_height,
)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=grid_savepoint.vct_a(),
vct_b=grid_savepoint.vct_b(),
)
vertical_params = v_grid.VerticalGrid(config=vertical_config, vct_a=grid_savepoint.vct_a())

metric_state = graupel.MetricStateIconGraupel(
ddqz_z_full=metrics_savepoint.ddqz_z_full(),
Expand Down
61 changes: 20 additions & 41 deletions model/common/src/icon4py/model/common/grid/vertical.py
Original file line number Diff line number Diff line change
Expand Up @@ -126,32 +126,25 @@ class VerticalGrid:
"""
Contains vertical physical parameters defined on the vertical grid derived from vertical grid configuration.

_vct_a and _vct_b: See docstring of get_vct_a_and_vct_b. Note that the height index starts from the model top.
_vct_a: See docstring of get_vct_a. Note that the height index starts from the model top.
_end_index_of_damping_layer: Height index above which Rayleigh damping of vertical wind is applied.
_start_index_for_moist_physics: Height index above which moist physics and advection of cloud and precipitation variables are turned off.
_end_index_of_flat_layer: Height index above which coordinate surfaces are flat.
"""

config: VerticalGridConfig
vct_a: dataclasses.InitVar[fa.KField[ta.wpfloat]]
vct_b: dataclasses.InitVar[fa.KField[ta.wpfloat]]
_vct_a: fa.KField[ta.wpfloat] = dataclasses.field(init=False)
_vct_b: fa.KField[ta.wpfloat] = dataclasses.field(init=False)
_end_index_of_damping_layer: Final[gtx.int32] = dataclasses.field(init=False)
_start_index_for_moist_physics: Final[gtx.int32] = dataclasses.field(init=False)
_end_index_of_flat_layer: Final[gtx.int32] = dataclasses.field(init=False)

def __post_init__(self, vct_a, vct_b):
def __post_init__(self, vct_a):
object.__setattr__(
self,
"_vct_a",
vct_a,
)
object.__setattr__(
self,
"_vct_b",
vct_b,
)
vct_a_array = self._vct_a.asnumpy()
object.__setattr__(
self,
Expand Down Expand Up @@ -248,10 +241,6 @@ def end_index_of_damping_layer(self) -> gtx.int32:
def vct_a(self) -> fa.KField:
return self._vct_a

@property
def vct_b(self) -> fa.KField:
return self._vct_b

def size(self, dim: gtx.Dimension) -> int:
assert dim.kind == gtx.DimensionKind.VERTICAL, "Only vertical dimensions are supported."
match dim:
Expand Down Expand Up @@ -293,11 +282,11 @@ def _determine_end_index_of_flat_layers(
)


def _read_vct_a_and_vct_b_from_file(
def _read_vct_a_from_file(
file_path: pathlib.Path, num_levels: int, allocator: gtx_typing.FieldBufferAllocationUtil
) -> tuple[fa.KField, fa.KField]:
"""
Read vct_a and vct_b from a file.
Read vct_a from a file.
The file format should be as follows (the same format used for icon):
# k vct_a(k) [Pa] vct_b(k) []
1 12000.0000000000 0.0000000000
Expand All @@ -311,19 +300,17 @@ def _read_vct_a_and_vct_b_from_file(
file_path: Path to the vertical grid file
num_levels: number of cell levels
allocator: GT4Py field allocator
Returns: one dimensional vct_a and vct_b arrays.
Returns: one dimensional vct_a arrays.
Comment on lines 289 to +303
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe explain here why vct_b is ignored (even though we keep the file format including vct_b for interoperability with fortran)

"""
num_levels_plus_one = num_levels + 1
vct_a = np.zeros(num_levels_plus_one, dtype=ta.wpfloat)
vct_b = np.zeros(num_levels_plus_one, dtype=ta.wpfloat)
try:
with file_path.open() as vertical_grid_file:
# skip the first line that contains titles
vertical_grid_file.readline()
for k in range(num_levels_plus_one):
grid_content = vertical_grid_file.readline().split()
vct_a[k] = ta.wpfloat(grid_content[1])
vct_b[k] = ta.wpfloat(grid_content[2])
except OSError as err:
raise FileNotFoundError(
f"Vertical coord table file {file_path} could not be read."
Expand All @@ -334,16 +321,14 @@ def _read_vct_a_and_vct_b_from_file(
) from err
except ValueError as err:
raise ValueError(f"data is not float at {k}-th line.") from err
return gtx.as_field((dims.KDim,), vct_a, allocator=allocator), gtx.as_field(
(dims.KDim,), vct_b, allocator=allocator
)
return gtx.as_field((dims.KDim,), vct_a, allocator=allocator)


def _compute_vct_a_and_vct_b( # noqa: PLR0912 [too-many-branches]
def _compute_vct_a( # noqa: PLR0912 [too-many-branches]
vertical_config: VerticalGridConfig, allocator: gtx_typing.FieldBufferAllocationUtil
) -> tuple[fa.KField, fa.KField]:
) -> fa.KField:
"""
Compute vct_a and vct_b.
Compute vct_a.

When the thickness of lowest level grid cells is larger than 0.01:
vct_a[k] = H (2/pi arccos((k/N)^s) )^d N = num_levels s = stretch_factor in vertical configuration
Comment on lines +331 to 334
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and add here a line, something like if needed vct_b can be computed as vct_b = np.exp(-vct_a / 5000.0) (so we don't have to look it up somewhere else)

Expand Down Expand Up @@ -379,8 +364,8 @@ def _compute_vct_a_and_vct_b( # noqa: PLR0912 [too-many-branches]

Args:
vertical_config: Vertical grid configuration
backend: GT4Py backend
Returns: one dimensional (dims.KDim) vct_a and vct_b gt4py fields.
allocator: GT4Py allocator
Returns: one dimensional (dims.KDim) vct_a gt4py fields.
"""
num_levels_plus_one = vertical_config.num_levels + 1
if vertical_config.lowest_layer_thickness > 0.01:
Expand Down Expand Up @@ -516,42 +501,36 @@ def _compute_vct_a_and_vct_b( # noqa: PLR0912 [too-many-branches]
)
/ ta.wpfloat(vertical_config.num_levels)
)
vct_b = np.exp(-vct_a / 5000.0)

if not np.allclose(vct_a[0], vertical_config.model_top_height):
log.warning(
f" Warning. vct_a[0], {vct_a[0]}, is not equal to model top height, {vertical_config.model_top_height}, of vertical configuration. Please consider changing the vertical setting."
)

return gtx.as_field((dims.KDim,), vct_a, allocator=allocator), gtx.as_field(
(dims.KDim,), vct_b, allocator=allocator
)
return gtx.as_field((dims.KDim,), vct_a, allocator=allocator)


def get_vct_a_and_vct_b(
def get_vct_a(
vertical_config: VerticalGridConfig, allocator: gtx_typing.FieldBufferAllocationUtil
) -> tuple[fa.KField, fa.KField]:
) -> fa.KField:
"""
get vct_a and vct_b.
get vct_a.
vct_a is an array that contains the height of grid interfaces (or half levels) from model surface to model top, before terrain-following coordinates are applied.
vct_b is an array that is used to initialize vertical wind speed above surface by a prescribed vertical profile when the surface vertical wind is given.
It is also used to modify the initial vertical wind speed above surface according to a prescribed vertical profile by linearly merging the surface vertical wind with the existing vertical wind.
See init_w and adjust_w in mo_nh_init_utils.f90.
Comment on lines 519 to 520
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
It is also used to modify the initial vertical wind speed above surface according to a prescribed vertical profile by linearly merging the surface vertical wind with the existing vertical wind.
See init_w and adjust_w in mo_nh_init_utils.f90.

these lines also refer to vct_b

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(if you think they're useful, move them to the other comment with the expression for vct_b)


When file_name is given in vertical_config, it will read both vct_a and vct_b from that file. Otherwise, they are analytically derived based on vertical configuration.
When file_name is given in vertical_config, it will read both vct_a from that file. Otherwise, they are analytically derived based on vertical configuration.

Args:
vertical_config: Vertical grid configuration
backend: GT4Py backend
Returns: one dimensional (dims.KDim) vct_a and vct_b gt4py fields.
allocator: GT4Py allocator
Returns: one dimensional (dims.KDim) vct_a gt4py fields.
"""

return (
_read_vct_a_and_vct_b_from_file(
vertical_config.file_path, vertical_config.num_levels, allocator
)
_read_vct_a_from_file(vertical_config.file_path, vertical_config.num_levels, allocator)
if vertical_config.file_path
else _compute_vct_a_and_vct_b(vertical_config, allocator)
else _compute_vct_a(vertical_config, allocator)
)


Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -267,11 +267,7 @@ def test_diagnostic_update_after_saturation_adjustement(
dtime = 2.0

vertical_config = v_grid.VerticalGridConfig(icon_grid.num_levels)
vertical_params = v_grid.VerticalGrid(
config=vertical_config,
vct_a=grid_savepoint.vct_a(),
vct_b=grid_savepoint.vct_b(),
)
vertical_params = v_grid.VerticalGrid(config=vertical_config, vct_a=grid_savepoint.vct_a())
virtual_temperature_tendency = data_alloc.zero_field(
icon_grid, dims.CellDim, dims.KDim, allocator=backend
)
Expand Down
4 changes: 1 addition & 3 deletions model/common/tests/common/fixtures.py
Original file line number Diff line number Diff line change
Expand Up @@ -158,9 +158,7 @@ def metrics_factory_from_savepoint(
stretch_factor=stretch_factor,
rayleigh_damping_height=damping_height,
)
vertical_grid = vertical.VerticalGrid(
vertical_config, grid_savepoint.vct_a(), grid_savepoint.vct_b()
)
vertical_grid = vertical.VerticalGrid(vertical_config, grid_savepoint.vct_a())
factory = metrics_factory.MetricsFieldsFactory(
grid=geometry_source.grid,
vertical_grid=vertical_grid,
Expand Down
Loading