Abstract This study presents a method for calculating all‐sky, broadband radiative fluxes for high‐resolution atmospheric model domains, such as limited‐area numerical weather prediction models and Cloud System‐Resolving Models, that have many columns (i.e. ≫ 10 5 ). The aim is to replace the Independent Column Approximation (ICA) which applies 1D radiative transfer models (RTMs) to all columns; but usually after skipping many dynamical time steps due to the ICA's computational burden. The proposed method, referred to as Partitioned Gauss–Legendre Quadrature (PGLQ), begins by partitioning a domain into sub‐domains. This study is concerned just with cloudy columns so partitioning was done according to cloud‐top altitudes. Then, for each sub‐domain, sort columns from smallest to largest cloud water path , assign each column a sequence number, identify n G GLQ points, and apply 1D RTMs to the n G columns with resulting flux profiles assigned to nearby columns in the sorted sequence. Last, reposition columns back to the full domain. Marine boundary‐layer and deep convective cloud fields were used to assess PGLQ. For ratios of RTM calculations between ICA and PGLQ of f RT = /( n G · ) −1 ∈ [10 3 , 10 4 ], MBEs and RMSEs for net short‐wave (SW) surface fluxes scale close to ( n G · ) −1 , RMSEs being typically ∼200× smaller than domain‐average . Structure function analyses of indicate that PGLQ errors resemble white noise. Long‐wave (LW) fluxes perform less well due to sorting on but seem likely to be acceptable, RMSEs being ∼5–10× smaller than domain‐average . All MBEs are typically <±0.5 W/m 2 . Estimates of radiative heating rate profiles are unbiased, and while noisy at small scales, they are portrayed well at scales larger than typical individual clouds. Overall speed‐up of PGLQ relative to ICA, for serial computation as consider here, is expected to be >300×; likely to be comparable to parallel computation. This should allow application of RTMs at every dynamical time step.