  This paper is concerned with the proof of existence and numerical
approximation of large-data global-in-time Young measure solutions to
initial-boundary-value problems for multidimensional nonlinear parabolic
systems of forward-backward type of the form $\partial_t u - \mbox{div}(a(Du))
+ Bu = F$, where $B \in \mathbb{R}^{m \times m}$, $Bv \cdot v \geq 0$ for all
$v \in \mathbb{R}^m$, $F$ is an $m$-component vector-function defined on a
bounded open Lipschitz domain $\Omega \subset \mathbb{R}^n$, and $a$ is a
locally Lipschitz mapping of the form $a(A)=K(A)A$, where $K\,:\, \mathbb{R}^{m
\times n} \rightarrow \mathbb{R}$. The function $a$ may have a nonstandard
growth rate, in the sense that it is permitted to have unequal lower and upper
growth rates. Furthermore, $a$ is not assumed to be monotone, nor is it assumed
to be the gradient of a potential. Problems of this type arise in mathematical
models of the atmospheric boundary layer and fall beyond the scope of monotone
operator theory. We develop a numerical algorithm for the approximate solution
of problems in this class, and we prove the convergence of the algorithm to a
Young measure solution of the system under consideration.
