Numerical simulations of wind turbine arrays are an essential tool for wind resource assessment, wind turbine siting, optimizing the performance of a wind farm, and analysis of loads on wind turbine blades under realistic operating conditions. Such simulations require accurate representation of atmospheric flows, wind shear and turbulence, and effects of operating turbines on such flows. However, numerical models most often used for simulation of flows in wind turbine arrays either do not include all the atmospheric effects, but critically depend on initial and boundary conditions provided by numerical weather prediction (NWP) models and coupling with them, or include unrealistic periodic boundary conditions thus effectively representing an infinite wind turbine array. Therefore, a single numerical model capable of both capturing large-scale atmospheric effects and resolving wind turbine interactions would be preferable. We therefore implemented a generalized actuator disc model in the Weather Research and Forecasting (WRF) NWP model. This actuator disc model is based on the NREL S809 open source blade design. We then demonstrated the use of WRF in large-eddy simulation (LES) mode to simulate flow through an array of wind turbines, a subset of an off-shore wind farm, under near-neutrally stratified conditions. Our simulations show that WRF-LES with a generalized actuator disk model is capable of capturing the characteristics of wind turbine wakes and their interactions with consecutive rows of wind turbines. The generalized actuator disk model implemented in WRF-LES is capable of representing the drag and torque effects distributed over the turbine rotor disk resolved with grid-cell size of several meters. We therefore conclude that WRF-LES with a generalized actuator disk model can be effectively used for detailed wind resource assessment, wind turbine siting, and blade load estimation.