Yes, this works. Here's a sketch: I'll work with $v=u^2$, which is also subharmonic. We then have that $$ v(z) \le \int_{\partial R} v(\zeta)\, d\omega_z(\zeta; R) , \quad\quad\quad\quad (1) $$ where $\omega$ denotes harmonic measure for a rectangle of the form $$ -q\eta\le x\le q\eta, \quad 1/4-q\le y\le 3/4+q , $$ say (in other words, the RHS of (1) solves the Dirichlet problem on $R$ with boundary values $v$).
We now simply average this formula over $1/2\le q \le 1$. Then the RHS becomes comparable to an average over an area, and thus is $\lesssim C$ by our assumption.
The rigorous treatment requires some simple estimates on harmonic measure as we vary $R$, but it is intuitively obvious that these will be comparable to arc length on the vertical pieces of the boundary, especially if we think of harmonic measure as hitting probability for Brownian motion.