Being a bit late for the party, here is nevertheless a small answer. In fact, there is a rather explicit way to compute adjoints of every differential operator (any order) between (smooth, compactly supported) sections of vector bundles. Of course, this is (in general) not the Hilbert space adjoint, here you need a bit more analysis ;)
The main idea is to use a symbol calculus based on a covariant derivative. I explain here the scalar version (trivial vector bundles) but the whole thing can be done in full generality as well. First, you choose a covariant derivative, say torsion-free and you choose a density on your manifold in order to have an integration measure. As you probably know, the symmetrized covariant derivative allows you to establish a $C^\infty(M)$-linear bijection between symbols, i.e. smooth functions on the cotangent bundle being polynomial in the fiber directions, and differential operators. Note that this is a real bijection, not just taking into account the leading symbol. Of course, the symbol depends on the chosen covariant derivative.
In a second step you compute once and for all the adjoint of a differential operator $D$ with symbol $f$ by zillions of integrations by parts. The funny thing is that there is a fairly simple way how the symbol of the adjoint looks like. You need two ingredients for that:
First, the covariant derivative allows you to define a horizontal lift which in turn determines a maximally indefinite pseudo Riemannian metric on the cotangent bundle (horizontal spaces are in bijection to tangent spaces at the base point, vertical spaces are in bijection to the cotangent space, thus there is a natural pairing). This metric has a Laplace operator (better: d'Alembert operator) $\Delta$ with which you can act on the symbol $f$. In the flat situation this is just
\begin{equation}
\Delta_{\mathrm{flat}} = \frac{\partial^2}{\partial q^i \partial p_i}
\end{equation}
for a Darboux chart on $T^*M$ induced by a chart on $M$. In general, there are a couple of Christoffel symbols needed to make this globally defined ;)
Second, the density $\mu$ of your integration might not be covariantly constant. In any case, it defines a one-form by
\begin{equation}
\alpha(X) = \frac{\nabla_X \mu}{\mu},
\end{equation}
which is now be used to cook up a new differential operator on $T^*M$. You can lift $\alpha$ vertically to a vector field $F(\alpha)$ on $T^*M$, completely canonical.
Having these two ingredients, the adjoint of $D^*$ has the following symbol
\begin{equation}
f^* = \exp\left(\frac{1}{2i}(\Delta + F(\alpha))\right) \overline{f}.
\end{equation}
The prefactor in the exponential depends a bit on your conventions concerning the assignment of symbol to operator. With this formula it is typically really just a computation to get adjoints of all kind of operators. In many cases, you can chose your density to be covariantly constant, so $\alpha = 0$.
You can find all this in much detail in my book on Poisson geometry, based on some old papers on quantization of cotangent bundles in the late 90s (together with Bordemann and Neumaier).
In the case of interesting bundles, the formula is essentially the same: you only have to choose covariant derivatives for the two vector bundles in question and modify the Laplace operator accordingly. Then you can proceed in the same way. This generalization is in a paper of mine with Bordemann, Neumaier and Pflaum.