We develop a residual-based variational multiscale (RBVMS) method based on isogeometric analysis for large-eddy simulation (LES) of wind-driven shear flow with Langmuir circulation (LC). Isogeometric analysis refers to our use of NURBS (Non-Uniform Rational B-splines) basis functions which have been proven to be highly accurate in LES of turbulent flows (Bazilevs, Y., 2007, Comput. Methods Appl. Mech. Eng., **197 **, pp. 173–201). LC consists of stream-wise vortices in the direction of the wind acting as a secondary flow structure to the primary, mean component of the flow driven by the wind. LC results from surface wave-current interaction and often occurs within the upper ocean mixed layer over deep water and in coastal shelf regions under wind speeds greater than 3 m s−1 . Our LES of wind-driven shallow water flow with LC is representative of a coastal shelf flow where LC extends to the bottom and interacts with the sea bed boundary layer. The governing LES equations are the Craik-Leivobich equations (Tejada-Martínez, A. E., and Grosch, C. E., 2007, J. Fluid Mech., **576 **, pp. 63–108; Gargett, A. E., 2004, Science, **306 **, pp. 1925–1928), consisting of the time-filtered Navier-Stokes equations. These equations possess the same structure as the Navier-Stokes equations with an extra vortex force term accounting for wave-current interaction giving rise to LC. The RBVMS method with quadratic NURBS is shown to possess good convergence characteristics in wind-driven flow with LC. Furthermore, the method yields LC structures in good agreement with those computed with the spectral method in (Thorpe, S. A., 2004, Annu. Rev. Fluids Mech., **36 **, pp. 584 55–79) and measured during field observations in (D’Alessio, S. J., , 1998, J. Phys. Oceanogr., **28 **, pp. 1624–1641; Kantha, L., and Clayson, C. A., 2004, Ocean Modelling, **6 **, pp. 101–124).