Large-eddy simulation (LES) is a well-established numerical technique, resolv- ing the most energetic turbulent fluctuations in the planetary boundary layer. By averaging these fluctuations, high-quality profiles of mean quantities and turbulence statistics can be obtained in experiments with well-defined initial and boundary conditions. Hence, LES data can be beneficial for assessment and optimisation of turbulence closure schemes. A database of 80 LES runs (DATABASE64) for neutral and stably stratified planetary boundary layers (PBLs) is applied in this study to optimize first-order turbulence closure (FOC). Approx- imations for the mixing length scale and stability correction functions have been made to minimise a relative root-mean-square error over the entire database. New stability functions have correct asymptotes describing regimes of strong and weak mixing found in theoretical approaches, atmospheric observations and LES. The correct asymptotes exclude the need for a critical Richardson number in the FOC formulation. Further, we analysed the FOC quality as functions of the integral PBL stability and the vertical model resolution. We show that the FOC is never perfect because the turbulence in the upper half of the PBL is not generated by the local vertical gradients. Accordingly, the parameterised and LES-based fluxes decor- relate in the upper PBL. With this imperfection in mind, we show that there is no systematic quality deterioration of the FOC in the strongly stable PBL provided that the vertical modelresolution is better than 10 levels within the PBL. In agreement with previous studies, we found that the quality improves slowly with the vertical resolution refinement, though it is generally wise not to overstretch the mesh in the lowest 500 m of the atmosphere where the observed, simulated and theoretically predicted stably stratified PBL is mostly located.