@exporter.export
@preprocess_and_wrap()
@check_units('[pressure]')
def weighted_continuous_average(pressure, *args, height=None, bottom=None, depth=None):
    r"""Calculate weighted-continuous mean of an arbitrary variable through a layer.

    Layer top and bottom specified in height or pressure.

    Formula based on that from [Holton2004]_ pg. 76 and the NCL function ``wgt_vertical_n``

    .. math::  WCA = \frac{\int_{p_s}^{p_b} A dp}{\int_{p_s}^{p_b} dp}

    where:

    * :math:`WCA` is the weighted continuous average of a variable.
    * :math:`p_b` is the bottom pressure level.
    * :math:`p_s` is the top pressure level.
    * :math:`A` is the variable whose weighted continuous average is being calculated.

    Parameters
    ----------
    pressure : `pint.Quantity`
        Atmospheric pressure profile.

    args : `pint.Quantity`
        Parameters for which the weighted-continuous mean is to be calculated.

    height : `pint.Quantity`, optional
        Heights from sounding. Standard atmosphere heights assumed (if needed)
        if no heights are given.

    bottom: `pint.Quantity`, optional
        The bottom of the layer in either the provided height coordinate
        or in pressure. Don't provide in meters AGL unless the provided
        height coordinate is meters AGL. Default is the first observation,
        assumed to be the surface.

    depth: `pint.Quantity`, optional
        Depth of the layer in meters or hPa.

    Returns
    -------
    list of `pint.Quantity`
        list of layer mean value for each profile in args.

    Notes
    -----
    Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
    Since this function returns scalar values when given a profile, this will return Pint
    Quantities even when given xarray DataArray profiles.

    """
    # Split pressure profile from other variables to average
    pres_prof, *others = get_layer(
        pressure, *args, height=height, bottom=bottom, depth=depth
    )

    return [trapezoid(var_prof, x=pres_prof) / (pres_prof[-1] - pres_prof[0])
            for var_prof in others]
