Add support for Generalized vertical height coordinate (levelType 150)

[marfel]

Perhaps you might consider supporting levelType 150, which is used by the state-of-the art models by the DWD, namely ICON and COSMO-D2. The following implementation allows reading the models on native levels, although one could probably do a more complete job by creating a corresponding HybridFactory or some such. I also did not write any tests, obviously ;P

def generalized_vertical_coordinates(section, metadata):
    """
    Translate the section 4 generalized vertical coordinates (level type 150).

    Updates the metadata in-place with the translations.

    Reference GRIB2 Code Table 4.5 and documentation of COSMO-D2 and ICON model
    by DWD.
    Parameter HHL must be read in to reconstruct the height grid; we
    refrain from giving misleading height and pressure coordinates. Model
    level refers to full levels (e.g. level centers). Half levels (e.g. level
    boundaries) consequently are assigned non-integer values.

    Args:

    * section:
        Dictionary of coded key/value pairs from section 4 of the message.

    * metadata:
        :class:`collections.OrderedDict` of metadata.

    """
    typeOfSecondFixedSurface = section['typeOfSecondFixedSurface']
    if typeOfSecondFixedSurface not in [150, 101, _TYPE_OF_FIXED_SURFACE_MISSING]:
        msg = 'Product definition section 4 contains unsupported type ' \
            'of second fixed surface [{}]'.format(typeOfSecondFixedSurface)
        raise TranslationError(msg)

    scaleFactor = section['scaleFactorOfFirstFixedSurface']
    if scaleFactor != 0:
        msg = 'Product definition section 4 contains invalid scale ' \
              'factor of first fixed surface [{}]'.format(scaleFactor)
        raise TranslationError(msg)

    # Create the model level number scalar coordinate.
    # In COSMO-D2 this references a layer boundary for TKE and W, and the
    # layer center of the layer below the boundary in other cases, hence
    # we average over the two level numbers given
    scaledValue = section['scaledValueOfFirstFixedSurface']

    if typeOfSecondFixedSurface in [101, _TYPE_OF_FIXED_SURFACE_MISSING]:
        # Message contains the level definition field HHL itself
        # or is defined on the half level
        # => use half level as a float value
        model_level = float(scaledValue) - 0.5
    else:
        # Message contains physical variable on level center
        model_level = scaledValue
    coord = DimCoord(model_level, standard_name='model_level_number',
                     attributes=dict(positive='down'))
    metadata['aux_coords_and_dims'].append((coord, None))


def vertical_coords(section, metadata):
    """
    Translate the vertical coordinates or hybrid vertical coordinates.

    Updates the metadata in-place with the translations.

    Reference GRIB2 Code Table 4.5.

    Args:

    * section:
        Dictionary of coded key/value pairs from section 4 of the message.

    * metadata:
        :class:`collections.OrderedDict` of metadata.

    """
    if section['NV'] > 0:
        # Generate hybrid vertical coordinates.
        if section['typeOfFirstFixedSurface'] == 150:
            generalized_vertical_coordinates(section, metadata)
        else:
            hybrid_factories(section, metadata)
...

One might also need to add some GRIB2 mappings, but I'm not sure. Just in case, here is what I currently define in addition to the standard mappings:

    # additions for ARPEGE
    G2Param(2, 0, 19, 11): CFName('turbulent_kinetic_energy', None, 'J kg-1'),
    G2Param(2, 0, 1, 16): CFName('snowmelt', None, 'kg m-2'),
    G2Param(2, 0, 1, 52): CFName('total_precipitation', None, 'kg m-2'),
    G2Param(2, 0, 2, 23): CFName('wind_speed_of_gust_u_comp', None, 'kg m-2'),
    G2Param(2, 0, 2, 24): CFName('wind_speed_of_gust_v_comp', None, 'kg m-2'),
    G2Param(2, 0, 1, 6): CFName('evaporation', None, 'kg m-2'),
    G2Param(2, 0, 2, 37): CFName('surface_downward_northward_stress', None, 'kg m-1 s-1'),
    G2Param(2, 0, 2, 38): CFName('surface_downward_eastward_stress', None, 'kg m-1 s-1'),
    G2Param(2, 0, 3, 18): CFName('planetary_boundary_layer_height', None, 'm'),
    G2Param(2, 0, 4, 196): CFName('clear_sky_downward_solar_flux', None, 'W m-2'),
    G2Param(2, 0, 5, 196): CFName('clear_sky_downward_thermal_flux', None, 'W m-2'),
    # additions for RACEVAR
    G2Param(2, 0, 1, 8): CFName('total_precipitation', None, 'kg m-2'),
    # additions for COSMO-REA6
    G2Param(2, 0, 1, 14): CFName('convective_snow_fall', None, 'kg m-2'),
    G2Param(2, 0, 1, 76): CFName('convective_rain_rate', None, 'kg m-2'),
    G2Param(2, 0, 1, 77): CFName('large_scale_rain_rate', None, 'kg m-2 s-1'),
    G2Param(2, 0, 6, 25): CFName('duration_sunshine', None, 'h'),
    G2Param(2, 2, 0, 5): CFName('water_runoff', None, 'kg m-2'),
    G2Param(2, 2, 3, 18): CFName('soil_temperature', None, 'K'),