Using adjoint optimization and three-dimensional Reynolds-averaged Navier Stokes (RANS) simulations, we developed a new gradient-based approach for optimally siting wind turbines within utility-scale wind plants. By solving the adjoint equations of the flow model, the gradients needed for optimization were found at a cost that was independent of the number of control variables, thereby permitting optimization of large wind plants with many turbine locations. Moreover, compared to the common approach of superimposing prescribed wake deficits onto linearized flow models, the computational efficiency of the adjoint approach allowed the use of higher-fidelity RANS flow models that can capture nonlinear turbulent flow physics within a wind plant. The RANS flow model was implemented in the Python finite element package FEniCS and the derivation of the adjoint equations was automated within the dolfin-adjoint framework. Layout optimization was demonstrated for complex wind roses, including a full annual energy production (AEP) layout optimization over 36 inflow directions and 5 windspeed bins.