As a simple alternative to the standard eddy-diffusivity closure, a nonlinear subgrid-scale (SGS) flux model is introduced and implemented in simulations of a neutral atmospheric boundary layer and a stable atmospheric boundary layer. The new model computes the structure of the SGS flux (relative magnitude of the vector components) based on the normalized gradient vector, which is derived from the Taylor expansion of the exact SGS flux. The SGS magnitude is computed as the product of a SGS velocity scale and a SGS scalar concentration scale, which are estimated based on the local-equilibrium hypothesis. To resolve the instability issue of the original gradient model and ensure numerical stability, we adopt a clipping procedure to avoid local negative SGS dissipation rate of the scalar variance. The model formulation, using constant coefficients, is assessed through a systematic comparison with well-established theoretical predictions and reference results of various flow statistics. Simulation results obtained with the use of this new model show good agreement with the reference results and an evident improvement over results obtained using traditional eddy-diffusivity models. For instance, the new model can deliver the expected surface-layer similarity scalar profile and power-law scaling of the power spectrum of scalar fluctuation. (C) 2013 American Institute of Physics. [http://dx.doi.org/10.1063/1.4774342]