 You want to use “element_boundary” here. “skeleton” does not give you what you want.
“skeleton” terms are integrals over facets, but there are multiple volume elements attached to a facet.
The term
a += h*(ppbar)*(qqbar)*dx(skeleton=True)
is not well defined. Which of the volume elements attached to a facet should provide “p” and “q”?
With “element_boundary”, it is well defined. Every facet is integrated over twice  once as part of the boundary of each volume element.
“skeleton” is usually used when you have jumps between volume elements, something like
a += h*(pp.Other())*(qq.Other())*dx(skeleton=True)
Here, it does not matter which of the volume elements of a facet provides “p” and which one “p.Other()”.
To make this clear, your term, using “skeleton”, would be something like this:
a += h * ( (ppbar)*(qqbar) + (p.Other()pbar)*(q.Other()qbar) ) * dx(skeleton=True)
For “skeleton” to work, though, you have to set “dgjumps=True” in the FESpace:
X = FESpace([V, Vbar, Vbar, Q, Qbar], dgjumps=True)
Forgetting this, and using “skeleton” integrals leads to the illegal dnums error you encountered.

See the replies to this post from last week as for why your problem is not well posed if we use Dirichlet conditions everywhere and how to fix this issue.

There is a typo:
(u, ubar1, ubar2, p, pbar) = X.TrialFunction()
(v, vbar1, vbar2, q, qbar) = X.TrialFunction()
should be
(u, ubar1, ubar2, p, pbar) = X.TrialFunction()
(v, vbar1, vbar2, q, qbar) = X.TestFunction()
You can do
a.Assemble()
print(a.mat)
and look at the matrix entries as a first sanity check. In this case, ALL entries are 0, which leads us to suspect something is off with trial/test functions.
 Your formulation is wrong
bh_dx =  p*div(v) + q*div(u)
bh_ds = (vvbar)*n*pbar  (uubar)*n*qbar,
should, I think, be
bh_dx = p*div(v)  q*div(u)
bh_ds = v*n*pbar + u*n*qbar
You should multiply the divergence criterium by 1, or define p as p such that you get a symmetric matrix. This is not strictly speaking necessary, but good practice.
The term (vvbar)npbar is, I believe, wrong. I think it should be a bnpbar
Best,
Lukas
Attachment: sto.py