This paper is concerned with a first-order fully discrete decoupled method for solving the three-dimensional (3D) primitive equations of the ocean with the Dirichlet boundary conditions on the side, where a decoupled semi-implicit scheme is used for the time discretization, and the P1(P1)−P1−P1(P1) finite element for velocity, pressure, and density is used for the spatial discretization of these equations. The H1−L2−H1 optimal error estimates for the numerical solution (uhn,phn,θhn) and the L2 optimal error estimate for (uhn,θhn) are established under the restriction of 0<h≤β1 and 0<τ≤β2 for some positive constants β1 and β2. Moreover, numerical investigations are provided to show that the first-order decoupled method is of almost unconditional convergence with accuracy O(h+τ) in the H1-norm and O(h2+τ) in the L2-norm for solving the 3D primitive equations of the ocean. Numerical results are given to verify the theoretical analysis.