Abstract In convective flows, vertical turbulent fluxes, covariances between vertical velocity and scalar thermodynamic variables, include contributions from local mixing and large-scale coherent motions, such as updrafts and downdrafts. The relative contribution of these motions to the covariance is important in turbulence parameterizations. However, the flux partition is challenging, especially in regions without convective cloud. A method to decompose the vertical flux based on the corresponding joint probability density function (JPD) is introduced. The JPD-based method partitions the full JPD into a joint Gaussian part and the complement, which represent the local mixing and the large-scale coherent motions, respectively. The coherent part can be further divided into updraft and downdraft parts based on the sign of vertical velocity. The flow decomposition is independent of water condensate (cloud) and can be applied in cloud-free convection, the subcloud layer, and stratiform cloud regions. The method is applied to large-eddy simulation model data of three boundary layers. The results are compared with traditional cloud and cloud-core decompositions and a decaying scalar conditional sampling method. The JPD-based method includes a single free parameter and sensitivity tests show weak dependence on the parameter values. The results of the JPD-based method are somewhat similar to the cloud-core and conditional sampling methods. However, differences in the relative magnitude of the flux decomposition terms suggest that an objective definition of the flow regions is subtle and diagnosed flow properties like updraft characteristics depend on the sampling method. Moreover, the flux decomposition depends on the thermodynamic variable and convection characteristics.