In this paper, we investigate the robust linear precoder design for three dimensional (3D) massive multi-input multi-output (MIMO) downlink with uniform planar array (UPA) and imperfect channel state information (CSI). In practical massive MIMO with UPAs, the number of antennas in each column or row is usually limited. The straightforward extension of the conventional DFT based beam domain channel model widely used in massive MIMO with uniform linear arrays (ULAs) can not apply. To overcome this issue, we establish a new beam domain channel model by using sampled steering vectors. Then, a novel method to obtain the beam domain channel power matrices and the instantaneous beam domain channel coefficients is proposed, and an a posteriori beam domain channel model which includes the channel aging and the spatial correlation is established. On the basis of the a posteriori channel model, we consider the robust precoder design with the expected weighted sum-rate maximization under a total power constraint. By viewing the power constraint as a Riemannian manifold, we transform the constrained optimization problem into an unconstrained optimization problem on the Riemannian manifold. Then, we derive an iterative algorithm to obtain the optimal precoders by setting the Riemannian gradient of the objective function to zero. Furthermore, we propose a low complexity robust precoder design by replacing the expected rates in the objective function with their upper bounds. Simulation results show that the proposed precoders can achieve significant performance gain than the widely used regularized zero forcing (RZF) precoder and signal to leakage noise ratio (SLNR) precoder.