I think that actually one can construct such an example using a little harmonic function theory.
First, notice that by RMT and by symmetry there is a conformal map that extends to a homeomorphism from the closed unit disc to the closed square (and which is actually conformal everywhere outside the vertices), $F:\mathbb D \to S, F(0)=0, F(\pm 1, \pm i)$ being the vertices of the square in counterclockwise order corresponding to the circle ordering (so if one fixes the image of one vertex say $F(1)=(-1,-1)$, the others are fixed eg $F(i)=(1,-1)$ etc).
(One can write a formula for this as $F(z)=c\int_0^z\frac{dz}{\sqrt{1-z^4}}$ but that is not needed)
Consider the finite non-negative measure on the unit circle given by $0, \log 2, \log 3, \log 4$ on the four open arcs $(1,i), (i-1), (-1,-i), (-i,1)$ and by zero in the four points or if you want the absolute continuous one given by $d\mu=qdt$ where $q$ takes the given values on the four open arcs and is unimportant what finite value we give it at the $4$ roots of unity of order $4$
Let $u_1(re^{i\theta})=\frac{1}{2\pi}\int_0^{2\pi}\frac{1-r^2}{1-2r\cos (\theta-t)+r^2}d\mu(t)$ the Poisson transform of $d\mu$ which is harmonic, bounded and positive in the open disc and satisfying $u_1(re^{it}) \to q(t)$ non-tangentially outside the four special points; actually note that $u_1(0)=\frac{1}{2\pi}\int_0^{2\pi}d\mu(t)=\frac {\log 2+\log 3+\log 4}{4} >2$.
So considering $g=u+iv$ holomorphic in the unit disc, $h(z)=e^{g(z)}$ almost satisfies the required properties on the four arcs since $|h(z)|=e^{u(z)}$ and $|h(0)|=e^{\frac{\log 24}{4}}=24^{1/4}>2$ (however $h$ is not continuous on the boundary, though it is "almost" so) and then clearly $f(w)=h(F^{-1}(w))=e^{g(F^{-1}(w))}$ almost satisfies the requires properties on the square and $|f(0)|>2$;
But now it is clear that $h_1(z)=(1-\epsilon)h(rz)$ for $r$ close enough to $1$ and $\epsilon>0$ small will do (will satisfy the boundness properties on the arcs and will be continuous, even holomorphic on the closed disc) and then $|h_1(0)|>2$ if $1-r$ hence $\epsilon$ are small enough, so taking $f_1(w)=h_1(F^{-1}(w))$ solves the problem.
Note that since $F, F^{-1}$ are not conformal at the vertices, $f_1$ is only continuous on the closed square (though holomorphic outside the vertices) despite that $h_1$ is holomorphic on the closed unit disc
Best Answer
Your observation that $|f(z)|=M$ is useful. From the equality we get, if $z=e^{it}$, $$ f(z)=M\,e^{ig(t)},\quad g\in C([0,2\,\pi],\mathbb{R}),\quad g(0)=g(2\,\pi)\pmod {2\,\pi} $$ and $$ \int_0^{2\pi}e^{ig(t)}e^{it}\,dt=2\,\pi\,e^{i\alpha} $$ for some $\alpha\in[0,2\,\pi)$. From this $$ \Re(e^{-i\alpha}\int_0^{2\pi}e^{ig(t)}e^{it}\,dt)=\int_0^{2\pi}\Re(e^{i(g(t)+t-\alpha)})\,dt=2\,\pi. $$ Deduce that $$ \Re(e^{i(g(t)+t-\alpha)})=1\quad \forall t $$ and that $g(t)+t-\alpha=0$ for all $t$. I hope I have not given too much detail.