In this paper, we study the joint cooperative beamforming and user association scheme for multi-cell multi-user full-dimension massive multiple-input multiple-out (MIMO) networks, where each user equipment is dynamically associated with one of the base stations. We propose to maximize the network capacity by jointly optimizing the beamforming vectors and a binary vector that represents the user association decisions. The binary vectors are further transformed into binary beam association factors (BAFs) as indicators of both the beamforming and the user association decisions. A three-step Gaussian belief propagation (GaBP) based distributed solver is proposed to solve the non-convex NP-hard problem considered in this paper. By relaxing the binary BAFs, we first transform the optimization problem into a linear programming. Subsequently, GaBP is used to efficiently obtain a feasible solution on BAFs in parallel, and two mapping algorithms are proposed to achieve different performance-complexity tradeoffs. Simulation results show that the proposed cooperative beamforming methods significantly outperform the benchmarks in the literature in terms of the network capacity with a low computation complexity.