Your question is a bit vague, but let me try the following statement, which might be the kind of answer you are looking for: If $M$ is a manifold and $S\to M$ is a vector bundle over $M$ endowed with a connection $\nabla$, then a 'total differential equation' for sections $s$ of $S$ is an equation of the form $\nabla_X s = 0$ for all vector fields $X$. Then the Frobenius integrability conditions for this system are simply $\Omega^\nabla=0$, where $\Omega^\nabla$ is the curvature of $\nabla$, i.e., it is the $2$-form on $M$ with values in $\mathrm{End}(S)$ that is defined by
$$
\Omega^\nabla(X,Y)s = \nabla_X\nabla_Ys - \nabla_Y\nabla_Xs - \nabla_{[X,Y]}s
$$
for all vector fields $X$ and $Y$ on $M$. So, yes, the Frobenius integrability conditions are 'covariant', i.e., they can be formulated independently of coordinates.
Now, one may wonder, "But suppose I start with some linear total differential system
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b \tag{1}
$$
in local coordinates. Does this apply?" The answer is 'yes': One regards the $u^a$ as the components of a section $s = (u^a)$ of a trivialized vector bundle $S$ and the matrices $\Gamma_i = \bigl(\Gamma^a_{bi}(x)\bigr)$ as sections of $\mathrm{End}(S)$ and one defines a connection $\nabla$ on this trivial bundle by the rule
$$
\nabla_{\frac{\partial}{\partial x^i}}s = \frac{\partial s}{\partial x^i} - \Gamma_i\ s.
$$
Then the curvature of $\nabla$ vanishes if and only if the Frobenius conditions are satisfied for the system (1).
If one has an inhomogeneous linear system
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{bi}(x) u^b + \Lambda^a_i(x), \tag{2}
$$
one can convert it to a homogeneous system by replacing $\Lambda^a_i(x)$ by $\Lambda^a_i(x)u^0$ in the above equation and adjoining the equations
$$
\frac{\partial u^0}{\partial x^i} = 0,
$$
making it a homogeneous linear system in one more variable. (Effectively, one allows the $a$ and $b$ indices (which run from, say, $1$ to $s$) to attain the value $0$, and sets $\Gamma^0_{bi}= \Gamma^0_{0i}=0$ and $\Gamma^a_{0i}=\Lambda^a_i$.) Again, the Frobenius conditions are equivalent to the vanishing of the curvature of the connection associated to the augmented system (now in $s{+}1$ unknowns rather than $s$ unknowns).
If one goes all the way to a fully nonlinear total differential system, say,
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{i}(x,u), \tag{3}
$$
then one can't (usually) reformulate this as the vanishing of the covariant derivatives of a section of something. However, there is still a 'covariant' coordinate-free interpretation of the Frobenius condition in terms of a tensor, but the tensor is defined on the total space $S$ (which has the $x$'s and $u$'s as local coordinates). The point is then, that, on $S$, one has a plane field $D$ (sometimes called a 'distribution') defined by the Pfaffian equations
$$
du^a - \Gamma^a_{i}(x,u)\ dx^i = 0
$$
(summation on $i$ intended). Associated to this plane field $D\subset TS$, there is a natural skew-symmetric linear operator
$$
\Phi: D\times D \to TS{\bigl/}D
$$
i.e., a section of the bundle $\bigl(TS{\bigl/}D\bigr)\otimes \Lambda^2(D^\ast)$ over $S$ that vanishes if and only if the system (3) satisfies the Frobenius conditions. It is this operator that is the 'covariant' interpretation of the Frobenius integrability conditions. Of course, this operator is defined via the Lie bracket. On sections $X$ and $Y$ of $D$, it is just
$$
\Phi(X,Y) \equiv [X,Y]\ \mathrm{mod}\ D,
$$
so one gets the usual statement that the Frobenius integrability conditions are equivalent to the condition that the sections of $D$ be closed under Lie bracket.
Finally, to get back to when one can interpret the Frobenius conditions as the vanishing of the curvature of a connection, there is a rather general formulation of this that goes back all the way to Lie: Suppose that the system (3) can be written in the form
$$
\frac{\partial u^a}{\partial x^i} = \Gamma^a_{i}(x,u)
= \gamma_{i\sigma}(x) F^{a\sigma}(u), \tag{4}
$$
for some functions $\gamma_{i\sigma}$ of $x$ and $F^{a\sigma}$ of $u$ (where summation over a new index $\sigma$ is assumed), and, further, assume that the vector fields
$$
U^\sigma = F^{a\sigma}(u)\ \frac{\partial\ \ }{\partial u^a}
$$
form a Lie algebra $L$, i.e., that there exist constants $c_\sigma^{\rho\tau}=-c_\sigma^{\tau\rho}$ such that $[U^\rho,U^\tau] = c_\sigma^{\rho\tau}U^\sigma$. Then the system (4) is what is known as a system of Lie type.
To avoid degenerate situations, one can assume, without loss of generality, that the $U^\sigma$ are a basis of $L$. In this case, one has that the Frobenius integrability condition for (4) is equivalent to the vanishing of the curvature of the $L$-valued connection $1$-form
$$
\omega = U^\sigma\ \gamma_{i\sigma}(x)\ dx^i.
$$
The use of this form of the equations belongs in the theory of symmetry analysis of differential equations, about which there is a very large literature.
Of course, recognizing when (3) can be written in the form (4) is not a trivial matter; there are methods for doing this involving Cartan's Method of Equivalence, but that would take a while to explain.
Ad 1: Yes, there is. The formula is
$$\nabla^*(X^\flat \otimes u) = - \nabla_X u -\mathrm{div}(X) \cdot u,$$
as can easily seen by local computation. Here, $X$ is a vector field and $X^\flat$ is the dual one form w.r.t. the metric.
Note that unless you have a scalar product on the bundle $E$, too, the dual operator will be an operator $\Gamma^\infty(T^*M \otimes E^*) \longrightarrow \Gamma^\infty(E^*)$. The above formula holds in both cases (depending on which case you are in, $u$ will be a section of $E$ or $E^*$, and the connection on the right will be either your given connection on $E$ or the corresponding dual connection on $E^*$ (this connection is always well-defined by the formula $(\nabla u)(s) = d(u(s)) - u(\nabla(s))$).
A similar formula holds even in the case that you don't have a metric on $M$, in which case the dual operator will send sections of $TM \otimes E^* \otimes |\Lambda|$ to sections of $E^* \otimes |\Lambda|$, where $|\Lambda|$ is the density bundle.
Ad 2: Let me first say that in my opinion, one uses the word "formal" only to indicate that one doesn't bother about any functional analytic meaning of the word adjoint whatsoever, and one does not automatically want to say that there has to be any setting where there is an "actual" adjoint.
However, to give a more constructive answer: Yes, this operator is always an adjoint operator in the functional analytic sense.
First, this is true for elliptic differential operators on manifolds. On compact manifolds, they have a unique closed extension, whose adjoint operator is given on smooth functions by the formal adjoint.
In general, consider for example the space $\mathscr{D}=\mathscr{D}(M, E)$ of test sections (which is not a Fréchet space unless $M$ is compact!). It can be embedded into its dual $\mathscr{D^\prime}(M, E^*)$. Now any differential operator $A$ is a continuous operator on $\mathscr{D}$, and its adjoint $A^\prime$ is an operator on $\mathscr{D}^\prime$. However, on the elements $\mathscr{D}\subset \mathscr{D}^\prime$ $A^\prime$ acts just as the formal adjoint of $A$. In this sense, the formal adjoint is also an "actual" adjoint in some sense.
However, this is quite formal and not at all deep. It follows at once from the definitions of distributions etc.
Best Answer
If I understand your notation correctly, then your question is a bit confused, because $g^N$ has to be a symmetric matrix, so that "$***$" = "$*$". The condition that $g^N$ is block diagonal does not have to hold; it says that the tangent vector of the last coordinate, $\partial/\partial u^{m+1}$, is perpendicular to the surface $M$. On the other hand, there always exist local coordinates with this property. If you take any local coordinates for $M$, you can evolve them for a short time with the normal surface flow. You can even get the condition $B = 1$ in a local chart.
Also, there certainly is another way to get the covariant derivative on $M$ and its Christoffel symbol. Namely, if you apply the covariant derivative $\nabla^N$ to a tangent vector field $v$ on $M$ in some tangent direction $w$, you get a vector field $\nabla^N_w(v)$ on $M$ that does not have to be tangent. You should then just project this derivative $\nabla^N_w(v)$ orthogonally onto the tangent bundle $TM$. The orthogonal projection is a useful tensor field $P$ defined on the tangent bundle $TN$ restricted to $M$, and you can write an explicit expression for the covariant derivative $\nabla^M$, or the Christoffel symbol or even the curvature tensor, in terms of $\nabla^N$ and this tensor field $P$. Actually, I am not entirely sure that this method is algebraically all that different, but it is at least conceptually different.