[ewquon]

Two fixes under column_fractional_z_index: which column it bisects, and the dispatch that reaches it. The second closes #508 and is what a conformal-conic parent needs in order to get its own column's vertical; the first is a pre-existing bug found while writing the second.

1. Pick the column by nearest node, not by truncation

column_fractional_z_index chose the column to bisect with Base.unsafe_trunc, which floors. Its ii/jj are 1-based node indices, so flooring takes the column west and/or south of the query whenever the fractional part is at least a half.

The sampled height is then off by however much that level's height surface tilts across one parent cell. Measured on an ERA5 0.25° parent over the Southern Great Plains (45×39×37, terrain 32–1518 m, 2016-09-10 12Z), queried from a 1/24° child at 12 heights above ground — 592,860 points, each interpolating the parent's own height field through the full stencil, so (sampled − requested) is the error this function contributes:

	floor	nearest
mean height error, all queries	−1.424 m	+0.153 m
rms height error, all queries	3.242 m	2.338 m
mean / rms, ≥ 500 m AGL	−1.915 / 3.016 m	−0.010 / 1.515 m
mean / rms, ≤ 250 m AGL	−0.442 / 3.651 m	+0.479 / 3.437 m

74.5 % of queries change column. What the fix buys is removing the bias: flooring always picks west and south, so its error carries a sign, and over a nested run that is a systematic displacement of the parent's height field rather than noise that averages out.

Across all 592,860 queries the largest temperature change is 0.151 K and the largest wind change 0.402 m/s, and nothing below 250 m AGL moves more than 7 m in height.

The residual rms does not go to zero and cannot. Its reference is the height solved on the bilinearly blended column, so what remains is the cost of picking one column at all: _interpolate blends four columns for the field value while column_fractional_z_index bisects a single one, which makes the height surface zeroth-order where the field on it is first-order. Going first-order would mean bisecting a blended profile, and that sits on the hot path for every boundary getbc, relaxation target and child initial condition — worth its own issue rather than this one.

A second path reaches the same defect, and is what the lat-lon regression test exercises: where the grid spacing is not a binary fraction, a query placed exactly on a node round-trips a few ulp below the integer and floors a whole column down.

2. Compose the conformal-conic horizontal with the per-column vertical

Closes #508.

PressureLevelVerticalDiscretization specializes _fractional_indices on the vertical; LambertConformalConicGrid specializes it on the horizontal, because its projected axes are not separable and cannot route through fractional_x_index/fractional_y_index. On a grid that is both, the Lambert method is more specific and wins outright — no ambiguity, no error — and the vertical falls back to the generic fractional_z_index, which reads the column-mean height profile.

Forcing the two apart on an 8×6×5 grid whose column i spans 100i to 100i + 1600 m, so that the domain-mean profile (450–2050 m) matches no column. k is the fractional level index _fractional_indices returns, on 5 levels spaced 400 m apart, so k = 1.5 sits halfway between a column's first and second level rather than a fraction of its depth:

column	that column spans	query height	k before	k after	exact k
i = 1	100–1700 m	300 m	0.625	1.500	1.500
i = 3	300–1900 m	1000 m	2.375	2.750	2.750
i = 8	800–2400 m	1000 m	2.375	1.500	1.500

The two 1000 m rows are the tell: before, columns 3 and 8 return the same index although their heights differ by 500 m, because both read the domain mean rather than their own column. Row 1 shows the second hazard — k = 0.625 is below level 1, underground, and nothing clamps it; a query above a column's top runs past k > Nz the same way. column_fractional_z_index clamps both ends.

This is not specific to conformal conic grids: any grid that specializes _fractional_indices jointly on the horizontal defeats a vertical specialization the same way. A horizontal_fractional_indices hook upstream would let the two compose rather than compete, and the joint method here is the local fix. It imports lcc_fractional_indices, which Oceananigans does not export; reimplementing it from the exported lcc_forward would duplicate the Center/Face node offsets.

Verification

On Oceananigans 0.114.0 and Julia 1.12.7:

• test/test_pressure_level_grid.jl passes 1466/1466 on CPU under --check-bounds=yes, including two new testsets. The lat-lon one fails without fix 1; the Lambert one fails without fix 2.
• No new method ambiguities. Pairwise Base.isambiguous over all 25 _fractional_indices methods finds 8, every one of them pre-existing and none involving the new method. The Union{Center, Face} restriction is load-bearing: dropping it adds eight more.
• Type-stable and allocation-free — @inferred concrete, @allocated 0 — and it runs inside a KernelAbstractions kernel on the CPU backend.
• Unaffected paths dispatch where they did before: plain LatitudeLongitudeGrid and plain LambertConformalConicGrid reach Oceananigans' own methods, and LatitudeLongitudeGrid + PressureLevelVerticalDiscretization reaches the existing method in this file.

Compatibility

No API change: no signatures or exports move, and the new _fractional_indices method is additive.

Results do change for every PressureLevelVerticalDiscretization parent, by the amounts in §1 — at most 0.151 K and 0.402 m/s on the grid measured there. Nested runs will not reproduce bit-for-bit across this commit. Flagging rather than deciding: whether that warrants a minor bump or ships as a fix is a maintainer call.

[codecov]

Codecov Report

❌ Patch coverage is 87.50000% with 1 line in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
...rc/Grids/pressure_level_vertical_discretization.jl	87.50%	1 Missing ⚠️

📢 Thoughts on this report? Let us know!

[glwagner]

Suggested change

	const LCCPressureLevelGrid =
	LambertConformalConicGrid{<:Any, <:Any, <:Any, <:Any, <:PressureLevelVerticalDiscretization}
	const LCCPressureLevelGrid = LambertConformalConicGrid{<:Any, <:Any, <:Any, <:Any, <:PressureLevelVerticalDiscretization}