Let $\gamma_R$ denote the positively oriented contour $[-R,R]$ followed by $Re^{it},0\le t\le \pi.$ If $z$ is in the upper half plane, then for $R>|z|$ we can apply Cauchy's formula to get
$$f(z) = \frac{1}{2\pi i}\int_{\gamma _R} \frac{f(w)}{w-z}\,dw$$
Now let $R\to \infty$ and use the assumed decay rate of $f$ to see
$$f(z) = \frac{1}{2\pi i}\int_\infty^\infty \frac{f(t)}{t-z}\,dt.$$
Claim: If $f$ is replaced by $\bar f$ in the last integral, the integral equals $0.$ Assuming the claim, we have
$$f(z) = \frac{1}{2\pi i}\int_\infty^\infty \frac{f(t)+\bar f (t)}{t-z}\,dt=\frac{1}{2\pi i}\int_\infty^\infty \frac{2\text {Re}f(t)}{t-z}\,dt=\frac{1}{\pi i}\int_\infty^\infty \frac{\text {Re}f(t)}{t-z}\,dt.$$
That's the desired result. I'll leave the proof of the claim for you to ponder for now.
No, the conjugate of a harmonic function that is continuous up to the boundary need not be continuous up to the boundary, or even bounded. This is related to the fact that the Hilbert transform does not preserve continuity (though it does preserve Hölder continuity of exponents $\alpha\in (0,1)$).
Here is an example. Let $F$ be a conformal map of the unit disk onto the domain bounded by the curves $y=0$ and $y=1/(1+x^2)$. Clearly, $F$ is not bounded: there are two points on the boundary of the unit disk which are sent into infinity by $F$. However, the imaginary part of $F$ is continuous: at the aforementioned points it approaches $0$. Thus, $\operatorname{Im}F$ is a harmonic function that is continuous in the closed unit disk, but whose conjugate is not bounded in the unit disk.
Another example is given by conformal map $G$ onto this rectangle with infinitely many slits:
Here the real part of $G$ is continuous up to the boundary, but the imaginary part is discontinuous, despite being bounded. (There is a point of unit circle where its cluster set is an interval of the size equal to the height of these slits.)
Behaviour of the convolution with the Poisson kernel near the discontinuous points
Typical example: consider the behavior of $\arg z$ in the upper halfplane as $z$ approaches $0$. Note that the boundary values of $\arg z$ have a jump discontinuity at $0$: they jump from $0$ to $\pi$. When $z$ approaches $0$ from within the domain, we can get any number between $0$ and $\pi$ as a limit. Or have no limit at all, if the curve of approach is zig-zagging.
This behavior is common. In fact, one can use $\arg z$ to prove this claim: subtracting an appropriate version of it from the boundary values removes the discontinuity, and then the behavior of function becomes clear: it consists of continuous part, plus a part that behaves like $\arg z$ above.
Best Answer
Consider $f_\alpha (z):=f(\alpha z)$ where $0<\alpha<1$. Now $f_\alpha $ is analytic on the unit circle.Use the Cauchy integral formula to obtain the following: $$ I_\alpha:=\frac1{2\pi} \int_0^{2\pi} f_\alpha(e^{i\theta}) P(r,\theta-t) d\theta =f(\alpha a)$$ As $\alpha \to 1$ , $f(\alpha a) \to f(a)$ and $I_\alpha \to I$ ( by the boundedness of $P$ and uniform convergence of $f_\alpha $ to $f$ on the unit circle)
So we get the Poisson formula.