Three-dimensional numerical simulation of equilibrium micromagnetic configurations existing in thin ferromagnetic films with surface anisotropy is carried out taking into account the strong demagnetization field acting on the film magnetization and the true micromagnetic boundary condition on the film surface. The numerical results are obtained in the simplest Neel approximation for surface anisotropy energy, a surface anisotropy constant K-s being a single phenomenological parameter. It is found that the spin canted state has the lowest total energy as compared to various multi-domain configurations in the intermediate range of thickness, L-z,(min) < L-z < L-z,L-max, if the magnitude of surface anisotropy constant Ks is below a certain critical value. For small thickness, L-z < L-z,L-min, the film is perpendicular magnetized, whereas for a thicker film, L-z > L-z,L-max, nearly uniform in-plane magnetization, or the vortex has been obtained depending on the film in-plane aspect ratio. On the other hand, different labyrinth domain structures with large in-plane magnetization have been calculated in a thick enough film, L-z > L-z,(max), with a sufficiently large surface anisotropy constant.