The comment I wrote above is perhaps too sketchy to be useful to the OP, so here is a somewhat expanded answer:

**Initial Remark:** If one doesn't require that the planar coordinate curves be orthogonal, then every surface has such coordinates locally, and one can construct them by simply slicing the surface with two generic $1$-parameter families of planes. Thus, from now on, I'll concentrate on the more interesting problem of understanding the surfaces that possess local orthogonal foliations by planar curves.

One can think of the problem of parametrizing such surfaces as that of finding the mappings $X:\mathbb{R}^2\to\mathbb{R}^3$ (where $u$ and $v$ are coordinates on the domain) that satisfy the three differential equations
$$
X_u\cdot X_v = 0,
\qquad X_u\wedge X_{uu}\wedge X_{uuu} = 0,
\qquad X_v\wedge X_{vv}\wedge X_{vvv} = 0.
$$
For simplicity, one should probably also require the *nondegeneracy conditions*
$$
X_u\wedge X_{uu} \not= 0,
\qquad X_v\wedge X_{vv} \not = 0.
$$
so that each coordinate curve lies in a unique $2$-plane in $\mathbb{R}^3$ and has nonvanishing curvature there.

Because this is three partial differential equations (one of order $1$ and two of order $3$), for the three unknown components of $X(u,v)$, this is at least formally determined, and so it should not be surprising that there are solutions.

However, it turns out that the symbol of this system is somewhat degenerate; for example, how one might write the system in Cauchy form is not immediate, and it's not clear what data should be provided for a well-posed Cauchy problem.

Now, a little calculation with the moving frame and exterior differential systems shows that Cauchy data are properly provided along a curve in the $uv$-plane on which both $du$ and $dv$ are nonzero. In fact, both $du=0$ and $dv=0$ are (not surprisingly) characteristics for the prolonged equation, and each of these has multiplicity $3$. It turns out, though, that, for each solution $X(u,v)$ there is a *seventh* characteristic foliation, this one depending on the particular solution $X(u,v)$ and being distinct from the $du=0$ and $dv=0$ directions. This system, properly prolonged, becomes involutive, and a solution is determined by specifying a certain set of $7$ essentially arbitrary functions along a noncharacteristic curve in the $uv$-plane.

Thus, the (local) solution surfaces that have such orthogonal planar coordinate curves depend on $5$ arbitrary functions of $1$ variable, where I am counting as equivalent all of the solutions $X\bigl(f(u),g(v)\bigl)$ where $f$ and $g$ are functions of a single variable, since the image surface is the same in these cases. *In particular, the generic surface in $\mathbb{R}^3$ does not, even locally, support a system of orthogonal planar coordinate curves.*

Finally, the OP asked about an algorithm for constructing these coordinates when they exist.
There is such an algorithm, and, when it is applied in a particular case (aside from a few degenerate cases that I'll describe below), the algorithm effectively either reduces the problem to constructing the integrals of a Frobenius system or returns an answer that there are no such coordinates. Here is a rough description of it:

First, if the given surface is a plane, then, obviously, any orthogonal coordinates on the surface will do, so set this case aside.

Second, if the surface is not a plane, but is ruled, then you could try using the ruling as one family of coordinate level sets, and you only need to check whether the orthogonal curve foliation consists of planar curves, which can be done by simple differentiation operations if the surface is explicitly given.

Third, if the surface is not ruled, then it is easy to show that neither family of coordinate curves of such a desired coordinate system can be asymptotic curves. In this case, we can think of the problem as specifying an orthogonal frame field (essentially a section of a circle bundle over the surface) whose legs are tangent to the two coordinate lines and are not pointed in any asymptotic directions.

Thinking of this as specifying an angle $\theta$ subject to some forbidden values, then the condition is that the two orthogonal curve families of the frame field be planar, which constitutes two second order PDE for the unknown function $\theta$. Because we are assuming that the legs are not asymptotic, one can show that the symbols of these two PDE have no common factor. Hence the prolongation of this system becomes a total system of PDE for the all of the third derivatives of $\theta$. (In fact, the system determines all but one of the second derivatives, so that the tableau of the prolongation vanishes at a point when the Pfaffian system has rank $4$.)

At this point, the standard methods for treating such total PDE systems reduce the problem to finding the integral surfaces of a $2$-plane field on an $6$-manifold. If one knows the surface explicitly (or merely knows the surface well enough to describe all of the relations on its second fundamental form and its covariant derivatives), then this is a differential algebra problem of finite type, and hence can be solved algorithmically.

For the generic surface, there will be no solutions, and the generic surface that does have a solution will have only one. Some special surfaces, such as the umbilic $2$-sphere, will have a finite dimensional family of solutions of positive dimension, but these surfaces will be quite rare. There are some obvious families of surfaces that have single solutions, such as the family of surfaces of revolution. However, I think it is not likely that there will be a simple characterization of the general surface that supports such a coordinate system.