Let $(\Sigma,g)$ be a compact Riemann surface with smooth boundary $\partial\Sigma$, $\Delta_g$ be the Laplace-Beltrami operator, and $h$ be a positive smooth function. Using a min-max scheme introduced by Djadli-Malchiodi (2006) and Djadli (2008), we prove that if $\Sigma$ is non-contractible, then for any $\rho\in(8k\pi,8(k+1)\pi)$ with $k\in\mathbb{N}^\ast$, the mean field equation $$\left\{\begin{array}{lll} \Delta_g u=\rho\frac{he^u}{\int_\Sigma he^udv_g}&{\rm in}&\Sigma\\[1.5ex] u=0&{\rm on}&\partial\Sigma \end{array}\right.$$ has a solution. This generalizes earlier existence results of Ding-Jost-Li-Wang (1999) and Chen-Lin (2003) in the Euclidean domain. Also we consider the corresponding Neumann boundary value problem. If $h$ is a positive smooth function, then for any $\rho\in(4k\pi,4(k+1)\pi)$ with $k\in\mathbb{N}^\ast$, the mean field equation $$\left\{\begin{array}{lll} \Delta_g u=\rho\left(\frac{he^u}{\int_\Sigma he^udv_g}-\frac{1}{|\Sigma|}\right)&{\rm in}&\Sigma\\[1.5ex] \partial u/\partial{\mathbf{v}}=0&{\rm on}&\partial\Sigma \end{array}\right.$$ has a solution, where $\mathbf{v}$ denotes the unit normal outward vector on $\partial\Sigma$. Note that in this case we do not require the surface to be non-contractible.