<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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()