NOAA NCEP Grib2 GFS & HRRR: levels, steps and duplicate variables!

It seems NCEP add several custom encodings to the WMO Grib2 standard that the CFGRIB library can't decode.

In the HRRR SubHourly product the Step is encoded as a range which breaks the CFGrib reader. When reading these variables with scan_grib you will get messages like:
2024-01-09T16:59:33.060Z MainProcess MainThread WARNING:grib2-to-zarr:Ignoring coordinate 'step' for varname 'vbdsf', raises: eccodes.WrongStepUnitError(Wrong units for step (step must be integer)) for variables: dswrf, vbdsf , tp , sdwe and unknown. Some of the grib messages really fail to decode - resulting in the unknowns. The step can be inferred from the runtime and the validtime of the model, but I think NCEP was trying to encode the duration of the average for the variables with stepType avg.

By comparing the results of using scan_grib and parsing the idx files provided by ncep, I was able to identify a few more edge cases. The table below shows some of the variables from gs://global-forecast-system/gfs.20231001/00/atmos/gfs.t00z.pgrb2.0p25.f006 which have duplicate variable name, step type, level type and level. Currently, the grib_tree method assumes these will be unique and silently takes the data from the last message in the file.

There are two types of duplicates I have found so far:

1. The GFS grib2 files include two accumulations for Convective Precipitation and Total Precipitation. One is the accumulation during the current model step and one is the total accumulation during the forecast run so far. With the step value parsed by CFGrib this is ambiguous for all model horizons (0 to 240 hour forecast files). With the idx file (gs://global-forecast-system/gfs.20231001/00/atmos/gfs.t00z.pgrb2.0p25.f006.idx) we can see a bit more metadata ACPCP:surface:0-6 hour acc fcst but for the first few timesteps of the model, even the idx values appear to be duplicates because the total is equal to the step accumulation.
2. There are several other variables that have level range such as 180-0 mb above ground and 0.44-1 sigma layer which decode as NaN with CFGrib (via kerchunk scan_grib). These result in additional duplicates which can confuse grib_tree (and anybody using it).

Fixing the actual decoding of the variables is hard. It may be possible by adding custom ecCodes definitions.

In the mean time, I want this issue to exist in the world for anyone also wondering what is going on.

Suggestions on improving the behavior of grib_tree in the mean time would be welcome. At present it is silently picking the last grib message and using the data (offset and length) for the given variable. This might be more than a little surprising for some users.

Start by reproducing the HRRR and GFS cases with scan_grib and grib_tree, using the linked product and index examples in the issue. Inspect how duplicate variable keys and unsupported level or step encodings are handled. Done should include an agreed behavior that no longer silently selects an ambiguous GRIB message, with coverage for the reported duplicates.