@sodoesaburningbus thanks for submitting this bug report. The error is in the static stability parameter that is the default value doesn't contain any units, even though the documentation states that the parameter should be a value with units, the default value violates that documentation!

@kgoebber Good catch and thanks @sodoesaburningbus for the report. I'll note that the reason this isn't currently erroring is because the unit-checking decorator isn't being given a dimensionality for static stability. I'll knock out a quick patch.

The current test test_q_vector_without_static_stability is actually specifically looking for these units. I think it was an oversight when we merged it 6 years ago and should have made the default correct on units. But that explains how this has passed as is.

I needed to check back with #661 and the textbooks to confirm this, but from what I recalled the unit difference was intended behavior, as it was accounting for differing definitions of Q-Vector (between Holton and Bluestein...Holton doesn't include the 1/sigma factor). Of course, over-6-year-ago me didn't clarify this in the docstring, so the mismatch in expectations is certainly understandable!

Though, @sodoesaburningbus, I did want to check: you mentioned expected units as kg/m2/s2. Our tests for the full version (accounting for static stability) have kg/m2/s3. Were you indeed expecting kg/m2/s2? If so, could you point us to a reference with that information, since our prior sources led (at least me) to believe it should be kg/m2/s3

I have #3690, but I do think we need to think hard about the default for the static stability parameter to q_vector. I don't love changing units, but I'm also not convinced that passing in a value of 1 with proper units is a good idea either. Unfortunately, changing the value starts to make this a bigger impact.

I now also took a look at Holton (4th ed.) and what is done there is keeping the static stability parameter on the operator side of the equation, which is different than what was done in the traditional QG-omega equation. Why he didn't divide that through in the same fashion to define Q is left unknown. In the Holton and Hakim (5th ed.), Q is derived in a substantially different manner (as is done in other parts of the text as well) from what had been the case through the first four editions. The equation is geared more toward a constant height as opposed to a constant pressure derivation.

Lackmann (2011) follows Bluestein, whereas Martin (2005) follows Holton (4th ed). So clear as mud. I think the leaving of the static stability parameter on the operator side is more done of convenience in writing the equation. In the end, static stability is one of the largest factors in modulating rising motion that is not either vorticity advection, temperature advection, or diabatic processes. As a result, this would be why I would advocate for using static stability with units in our function as opposed to not using it. Additionally, since we publish the Bluestein equation, which includes the static stability parameter, it seems logical that we would produce the calculation with the appropriate units for that particular equation.

It might also be worth checking to see how GEMPAK handles (or doesn't handle) the static stability parameter, if we are looking for feature parity. Or perhaps an optional argument to choose which flavor of the Q-vector to compute.

I've looked into this in more detail. First, GEMPAK has two Q-vector functions, QVEC, which doesn't include the static stability parameter, and QVCL, which does. Second, in the original Hoskins et al. (1978) paper that introduced Q-vectors (which uses a similar vertical coordinate as Holton/Hakim 5th ed), the static stability parameter ($N^2$ in this coordinate) is kept on the left-hand side (their Eq. 7). Based on my quick skim, Hoskins et al. don't discuss how static stability can modulate vertical motions at all.

Thus, the extant literature and software out there seems about equally divided. I think we should leave it up to the user. If they call q_vector and don't supply a value for static_stability, that means they want the version of Q without division by $\sigma$, with units to match. If they call q_vector and do supply a value of static_stability, that means they want the version of Q that does divide by $\sigma$, with units to match. I haven't actually tested the code yet, but I think this is what currently happens. If this is correct, then my vote would be "not a bug".

I would agree with @dopplershift that just tacking on units to the default value of 1 is not a good solution, since a value of 1 with SI units would be a pretty implausible value for $\sigma$ in the real world.

Regarding this comment:

The error is in the static stability parameter that is the default value doesn't contain any units, even though the documentation states that the parameter should be a value with units, the default value violates that documentation!

The documentation I see says "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." This seems to match the current implementation as far as I can tell.

@jthielen and I chatted during our MetPy call and I think we're in agreement with the last comment from @sgdecker and there isn't really a bug here. The ways I want to address this:

• Add missing unit-checking for the existing static_stability parameter (won't impact default at all)
• Leave the default of 1 (no units) to implement the equivalent GEMPAK QVEC and the matching textbook formulations
• Better document and cite the different formulations
• Better note the availability of static_stability as an alternative to the default (and match QVCL)
• Explicitly pass a representative value of static stability in our existing QVector.py example

@kgoebber @sodoesaburningbus Does that seem reasonable? Anything else to add to make this more clear for the future?