I may have to enter this as a sketch and fill in details later, but I thought that I'd go ahead and get the main ideas out there.
The first thing to notice is that the given problem is equivalent to the problem of solving $\rho(f)=g$ where $\rho:\mathrm{Imm}(T^2,\mathbb{R}^3)\to\Omega^2(T^2,\mathbb{R}^3)$ is the differential operator
$$
\rho(f) = f_x\times f_y\ \ dx{\wedge}dy
$$
taking values in (nonvanishing) $\mathbb{R}^3$-valued $2$-forms on $T^2$. The advantage of considering $\rho$ instead of $R$ is that $\rho$ is well-defined, independent of coordinates on $T^2$.
The second thing to notice is that the inner product on $\mathbb{R}^3$ is something of a red herring. The cross product is really just the wedge product followed by the orientation-and-metric induced canonical isomorphism between the second exterior power of $\mathbb{R}^3$ and $\mathbb{R}^3$ itself.
Let's undo all this, and consider, instead, a $3$-dimensional vector space $V$ and a surface $\Sigma$ and define an operator $\rho:\mathrm{Imm}(\Sigma,V)\to\Omega^2\bigl(T^2,\Lambda^2(V)\bigr)$ that, in local coordinates $(x,y)$ on $\Sigma$, has the expression
$$
\rho(f) = f_x{\wedge}f_y\ \ dx{\wedge}dy.
$$
Now we want to ask the question "Given a nonvanishing $\Lambda^2(V)$-valued $2$-form $g$ on a surface $\Sigma$, when can it be written in the form $g = \rho(f)$ for some immersion $f:\Sigma\to V$ and in how many ways?" This question is interesting both locally and globally. Locally, it is three (nonlinear) first-order equations for three unknowns, but, as we shall see, it is highly degenerate, even when $g$ is 'generic', so that, for example, there is no hope of writing it in Cauchy form.
One distinct advantage of formulating the equation in this way is that its full equivariance is manifest. Not only does $\rho$ not depend on a choice of local coordinates on the surface, if $L(v) = Av + b$ is any invertible affine transformation of $V$, then one has $\rho(L\circ f) = \Lambda^2(A)\bigl(\rho(f)\bigr)$. In particular, replacing $g$ by $\Lambda^2(A)(g)$ does not change the problem, which shows that one can use the geometry of linear transformations of $V$ to uncover invariants of $g$ that will determine properties of the solution.
For example, in the (very degenerate) case in which $g$ is a $2$-form that takes values in a line in $\Lambda^2(V)$, one can write $g = e_1{\wedge}e_2\ \Phi$ where $e_1$ and $e_2$ are linearly independent and $\Phi$ is (nonvanishing) $2$-form on $\Sigma$. Then the only possible solutions $f$ to $\rho(f)=g$ are to have $f$ immerse $\Sigma$ into a $2$-plane $P$ parallel to the subspace spanned by $e_1$ and $e_2$ in such a way that it pulls back the volume form on $P$ dual to $e_1{\wedge}e_2$ to be $\Phi$. This is now essentially one equation for two unknowns and is always locally solvable, though there may well be no global solution. This is a trivial special case, though, and doesn't give you any sense of the general case.
On this first pass, I will avoid discussing the singular and intermediate cases and go directly to the 'generic' case, the case in which the projectivization of $g$, i.e., $[g]:\Sigma\to \mathbb{P}\bigl(\Lambda^2(V)\bigr)\simeq\mathbb{RP}^2$ is an immersion. In this case, the problem reduces to either an elliptic or a hyperbolic problem, depending on the geometry of the mapping $g$.
Example: Before considering the general case, it's worth pausing to consider a simple example to illustrate where the analysis will take us: Let $\Sigma$ be the $xy$-plane and consider the two maps $f_\pm:\Sigma\to\mathbb{R}^3$ defined by
$$
f_\pm = \bigl(x,y,\tfrac12(x^2\pm y^2)\bigr).
$$
These have $g_\pm = \rho(f_\pm) = (1,0,x){\wedge}(0,1,\pm y)\ dx{\wedge}dy$. We want to find the immersions $f:\Sigma\to\mathbb{R}^3$ that satisfy $\rho(f)=g_\pm$. Well, such an $f$ will have to satisfy
$$
df = (1,0,x)\ \xi + (0,1,\pm y)\ \eta = (\xi,\eta,\ x\xi{\pm}y\eta\ )
$$
for some $1$-forms $\xi$ and $\eta$ that satisfy $\xi\wedge\eta=dx\wedge dy$ and are such that $d(df)=0$. Now, this latter equation implies $d\xi=d\eta=d(x\xi{\pm}y\eta)=0$, and so, assuming that a local solution is defined on a simply-connected open set $U\subset\Sigma$, one can then write $\xi = dp$ and $\eta = \pm dq$ for some functions $p$ and $q$ on $U$. The remaining equations on $p$ and $q$ then become $0=d(x\xi{\pm}y\eta)=dx\wedge dp + dy\wedge dq$ and $dp\wedge dq = \pm dx\wedge dy$. This is two PDE on $p$ and $q$ that can be written as a single Monge-Ampère equation by setting $p = u_x$ and $q = u_y$ (which works because $d(p\ dx + q\ dy) = 0$) and then the remaining equation becomes $u_{xx}u_{yy}-u_{xy}^2 = \pm 1$, the classic Monge-Ampère equation. This is, of course, elliptic for $g_+$ and hyperbolic for $g_-$. (We will see how this is reflected in the general analysis below.) The general solution to $\rho(f)=g_\pm$ is then
$$
f = (u_x,\ \pm u_y,\ xu_x+yu_y-u )
$$
where $u$ satisfies $u_{xx}u_{yy}-u_{xy}^2 = \pm 1$.
Now, I'm going to do the general analysis in terms of the moving frame and exterior differential systems, simply because that's the way I understand it and am most comfortable computing. Of course, one can avoid this, but I'm too lazy (and short on time) to do that now.
So fix a nondegenerate $g\in \Omega^2\bigl(\Sigma,\Lambda^2(V)\bigr)$. I'm going to define a bundle $P_g\to\Sigma$ whose elements will consist of the quadruples $(x;e_1,e_2,e_3)\in\Sigma\times V\times V\times V$ such that the $e_i$ are a basis of $V$ and, moreover, $e_1\wedge g(x) = e_2\wedge g(x) = 0$. On $P_g$, there exists a nonvanishing $2$-form $\Phi$ such that $g = e_1{\wedge}e_2\ \Phi$. This $\Phi$ is semi-basic for the projection to $\Sigma$ and hence is a multiple of any nonvanishing $2$-form on $\Sigma$ pulled up to $P_g$. Note that $P_g$ is a right $G_1$-bundle over $\Sigma$, where $G_1\subset\mathrm{GL}(3,\mathbb{R})$ is the subgroup that consists of matrices $A\in \mathrm{GL}(3,\mathbb{R})$ that are $(2,1)$-block upper triangular, i.e., $A^3_1=A^3_2=0$.
Now, one has the usual structure equations $de_a = e_b\ \eta^b_a$, which satisfy $d\eta^a_b = - \eta^a_c\wedge\eta^c_b$. (The summation convention always in force. I'll use the index ranges $1\le a,b,c \le 3$ and $1\le i,j,k \le 2$.) Since $g$ is a vector-valued $2$-form on a surface, one has $dg=0$, and expanding this using the structure equations yields
$$
0 = dg = e_1{\wedge}e_2\ \bigl(d\Phi + (\eta^1_1+\eta^2_2)\wedge\Phi\bigr)
+ e_3{\wedge}e_2\ \eta^3_1\wedge\Phi + e_1{\wedge}e_3\ \eta^3_2\wedge\Phi,
$$
so $\eta^3_1\wedge\Phi=\eta^3_2\wedge\Phi=0$. Now, the hypothesis that $[g]:\Sigma\to\mathbb{RP}^2$ be an immersion is equivalent to $\eta^3_1{\wedge}\eta^3_2$ being nonvanishing, Consequently, $\Phi = \lambda\ \eta^3_1{\wedge}\eta^3_2$ for some nonvanishing function $\lambda$ on $P_g$. One easily computes that, for $A\in G_1$,
$$
R_A^*\lambda = \left(\frac{A^3_3}{A^1_1A^2_2-A^2_1A^1_2}\right)^2\lambda,
$$
so the sign of lambda is constant on the fibers of $P_g$ and the equation $\lambda=\pm1$ defines a subset $P_g'\subset P_g$ that will be a right $G_2$-bundle over $\Sigma$, where $G_2\subset G_1$ is the subgroup consisting of those $A\in G_1$ that satisfy $A^3_3 = \pm(A^1_1A^2_2-A^2_1A^1_2)$.
I'll say that $g$ is of elliptic (resp. hyperbolic) type if $\lambda = +1$ (resp. $\lambda=-1$) on $P_g'$. (The reasons for these designations will become apparent below.)
Since $\Phi = \pm \eta^3_1{\wedge}\eta^3_2$ on $P_g'$, the structure equations plus the identity $0=d\Phi + (\eta^1_1+\eta^2_2)\wedge\Phi$ now show that $(\eta^3_3-\eta^1_1-\eta^2_2)\wedge\eta^3_1{\wedge}\eta^3_2=0$, so $\eta^3_3=\eta^1_1+\eta^2_2+b^1\eta^3_1+b^2\eta^3_2$ for some functions $b^1$ and $b^2$ on $P_g'$. The structure equations now show that the equations $b^1=b^2=0$ define a $G_3$-subbundle $P_g''\subset P_g'$, where $G_3\subset G_2$ is the subgroup consisting of those $A\in G_2$ with $A^1_3=A^2_3=0$. (This is the last reduction I will need to do to complete the calculation.)
Our goal now is to find $1$-forms $\eta^1$ and $\eta^2$ on $P_g''$ satisfying the equations
$$
d(e_1\ \eta^1 + e_2\ \eta^2) = 0\qquad\text{and}\qquad
e_1{\wedge}e_2\ \eta^1{\wedge}\eta^2 = g = \pm e_1{\wedge}e_2\ \eta^3_1{\wedge}\eta^3_2\ .
$$
If we can do this, then, at least locally, $e_1\ \eta^1 + e_2\ \eta^2 = df$ for some $V$-valued function $f:\Sigma\to V$ such that $\rho(f) = g$, and this will solve our problem.
Of course, the second equation is just $\eta^1\wedge\eta^2 = \pm\ \eta^3_1{\wedge}\eta^3_2$, while using the structure equations to expand the first equation yields
$$
0
= e_1 (d\eta^1 + \eta^1_1{\wedge}\eta^1+ \eta^1_2{\wedge}\eta^2)
+e_2 (d\eta^2 + \eta^2_1{\wedge}\eta^1+ \eta^2_2{\wedge}\eta^2)
+e_3 (\eta^3_1{\wedge}\eta^1 + \eta^3_2{\wedge}\eta^2)
$$
Looking at the $e_3$-coefficient in this equation and applying Cartan's Lemma, we see that there must exist $h^{ij}=h^{ji}$ such that $\eta^i = h^{ij}\ \eta^3_j$, which, substituted into the second equation, says that $\det(h) = \pm 1$.
Now, let $S_\pm$ denote the surface in the $3$-dimensional space of symmetric $2$-by-$2$ matrices that consists of the matrices with determinant $\pm 1$. ($S_+$ is a hyperboloid of 2 sheets, while $S_-$ is a hyperboloid of 1 sheet.) Let $h:P_g''\times S_\pm\to S_\pm$ be the projection onto the second factor. Set $\eta^i = h^{ij}\eta^3_j$ as $1$-forms on $P_g''\times S_\pm$ and consider the pair of $2$-forms
$$
\Upsilon^i = d\eta^i + \eta^i_1{\wedge}\eta^1+ \eta^i_2{\wedge}\eta^2.
$$
The formulae derived so far show that $\Upsilon^i = \pi^{ij}\wedge\eta^3_j$, where $\pi^{ij}=\pi^{ji}$ are $1$-forms defined by the relations
$$
\begin{pmatrix}\pi^{11}&\pi^{12}\\\\
\pi^{21}&\pi^{22}\end{pmatrix}
= dh + \begin{pmatrix}\eta^1_1&\eta^1_2\\\\
\eta^2_1&\eta^2_2\end{pmatrix} h
+ h\begin{pmatrix}\eta^1_1-\eta^3_3&\eta^2_1\\\\
\eta^1_2&\eta^2_2-\eta^3_3\end{pmatrix}.
$$
Moreover, the equation $\det(h)=\pm1$ and the normalizations made so far imply that $\mathrm{tr}(h^{-1}\pi)=0$.
These equations show that, on the quotient manifold $Q = (P_g''\times S_\pm)/G_3$, which is a bundle over $\Sigma$ with $2$-dimensional fibers, there is an ideal $\mathcal{I}$ that is locally generated by a pair of $2$-forms that pull back up to $P_g''\times S_\pm$ to be independent linear combinations of the $\Upsilon^i$, and that the integral surfaces of this ideal that are transverse to the fibers of $Q\to\Sigma$ represent solutions to our problem. When $\det(h)=1$, the above relations show that this is an elliptic ideal, while, when $\det(h)=-1$, this ideal is hyperbolic. In either case, the ideal is involutive, with Cartan characters $s_1=2$, $s_2=0$.
The Cartan-Kähler Theorem now applies, at least when $g$ is real-analytic, to show that there are local integral manifolds transverse to the fibers of $Q\to\Sigma$. In the hyperbolic case, a smooth version of the Cartan-Kähler theorem guarantees local solvability while, in the elliptic case, there are existence theorems (basically coming from the theory of pseudo-holomorphic curves) that guarantee existence of local solutions. Thus, the equations are locally solvable if $g$ is smooth.
Anyway, that's a sketch of an analysis of the nondegenerate case. The point is that the original determined system of three equations for three unknowns is not formally integrable, but, in the generic case, it can be prolonged to either an elliptic or a hyperbolic system that is involutive and hence formally integrable.
One should still do the case when the rank of the differential of $g$ is $1$, and one should look for a way to reinterpret the system $\mathcal{I}$ as a single scalar equation of Monge-Ampere type, to which more standard PDE methods would apply. There may even be a way to linearize this equation as either an elliptic or hyperbolic single scalar equation, but that would require more thought than I have had time to put into it so far.
