Define $F:\mathbb{R}_{x,y,z>0} \to \mathbb{R}$ by $F(x,y,z) = xy + xz + yz.$ Then $dF = (y+z, x+z, x+y)^T$ so $F$ has no critical values, and $1$ is a regular value of $F.$ The preimage of a regular value is a regular surface, so $S= F^{-1}(1)$ is a regular surface.
Then for points on the surface $xy+xz+yz=1$, we have $z = \frac{1-xy}{x+y}$ so we have to maximize $g(x,y)=\frac{ xy (1-xy)}{x+y}$ over $\mathbb{R_{x,y>0}}.$
We compute that $g_x = -\frac{y}{(x+y)^2} (x^2+2xy-1)$ and $g_y = -\frac{x}{(x+y)^2} (y^2+2xy-1)$ so $g_x=0$ when $y= \frac{1-x^2}{2x}$ and $g_y=0$ when $x= \frac{1-y^2}{2y}.$ When both these hold we have $1-x^2=2xy=1-y^2$ so $x=y$, and $x^2+2xy-1=0$ leads to $x=y=\frac{1}{\sqrt{3}}$ so the maximum value of $g$ is $\frac{1}{3\sqrt{3}}.$