def q_vector(u, v, temperature, pressure, dx, dy, static_stability=1):
    r"""Calculate Q-vector at a given pressure level using the u, v winds and temperature.

    .. math:: \vec{Q} = (Q_1, Q_2)
                      =  - \frac{R}{\sigma p}\left(
                               \frac{\partial \vec{v}_g}{\partial x} \cdot \nabla_p T,
                               \frac{\partial \vec{v}_g}{\partial y} \cdot \nabla_p T
                           \right)

    This formula follows equation 5.7.55 from [Bluestein1992]_, and can be used with the
    the below form of the quasigeostrophic omega equation to assess vertical motion
    ([Bluestein1992]_ equation 5.7.54):

    .. math:: \left( \nabla_p^2 + \frac{f_0^2}{\sigma} \frac{\partial^2}{\partial p^2}
                  \right) \omega =
              - 2 \nabla_p \cdot \vec{Q} -
                  \frac{R}{\sigma p} \beta \frac{\partial T}{\partial x}.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind (geostrophic in QG-theory)
    v : (M, N) ndarray
        y component of the wind (geostrophic in QG-theory)
    temperature : (M, N) ndarray
        Array of temperature at pressure level
    pressure : `pint.Quantity`
        Pressure at level
    dx : float or ndarray
        The grid spacing(s) in the x-direction. If an array, there should be one item less than
        the size of `u` along the applicable axis.
    dy : float or ndarray
        The grid spacing(s) in the y-direction. If an array, there should be one item less than
        the size of `u` along the applicable axis.
    static_stability : `pint.Quantity`, optional
        The static stability at the pressure level. Defaults to 1 if not given to calculate
        the Q-vector without factoring in static stability.

    Returns
    -------
    tuple of (M, N) ndarrays
        The components of the Q-vector in the u- and v-directions respectively

    See Also
    --------
    static_stability

    Notes
    -----
    If inputs have more than two dimensions, they are assumed to have either leading dimensions
    of (x, y) or trailing dimensions of (y, x), depending on the value of ``dim_order``.

    """
    dudy, dudx = gradient(u, deltas=(dy, dx), axes=(-2, -1))
    dvdy, dvdx = gradient(v, deltas=(dy, dx), axes=(-2, -1))
    dtempdy, dtempdx = gradient(temperature, deltas=(dy, dx), axes=(-2, -1))

    q1 = -mpconsts.Rd / (pressure * static_stability) * (dudx * dtempdx + dvdx * dtempdy)
    q2 = -mpconsts.Rd / (pressure * static_stability) * (dudy * dtempdx + dvdy * dtempdy)

    return q1.to_base_units(), q2.to_base_units()
r"""Calculate Q-vector at a given pressure level using the u, v winds and temperature.

    .. math:: \vec{Q} = (Q_1, Q_2)
                      =  - \frac{R}{\sigma p}\left(
                               \frac{\partial \vec{v}_g}{\partial x} \cdot \nabla_p T,
                               \frac{\partial \vec{v}_g}{\partial y} \cdot \nabla_p T
                           \right)

    This formula follows equation 5.7.55 from [Bluestein1992]_, and can be used with the
    the below form of the quasigeostrophic omega equation to assess vertical motion
    ([Bluestein1992]_ equation 5.7.54):

    .. math:: \left( \nabla_p^2 + \frac{f_0^2}{\sigma} \frac{\partial^2}{\partial p^2}
                  \right) \omega =
              - 2 \nabla_p \cdot \vec{Q} -
                  \frac{R}{\sigma p} \beta \frac{\partial T}{\partial x}.

    Parameters
    ----------
    u : (M, N) ndarray
        x component of the wind (geostrophic in QG-theory)
    v : (M, N) ndarray
        y component of the wind (geostrophic in QG-theory)
    temperature : (M, N) ndarray
        Array of temperature at pressure level
    pressure : `pint.Quantity`
        Pressure at level
    dx : float or ndarray
        The grid spacing(s) in the x-direction. If an array, there should be one item less than
        the size of `u` along the applicable axis.
    dy : float or ndarray
        The grid spacing(s) in the y-direction. If an array, there should be one item less than
        the size of `u` along the applicable axis.
    static_stability : `pint.Quantity`, optional
        The static stability at the pressure level. Defaults to 1 if not given to calculate
        the Q-vector without factoring in static stability.

    Returns
    -------
    tuple of (M, N) ndarrays
        The components of the Q-vector in the u- and v-directions respectively

    See Also
    --------
    static_stability

    Notes
    -----
    If inputs have more than two dimensions, they are assumed to have either leading dimensions
    of (x, y) or trailing dimensions of (y, x), depending on the value of ``dim_order``.