Cleo: the fundamental design of a new superdroplet model for high computational performance (v0.39.0)

Cleo is a novel implementation of the Super-Droplet Method (SDM) aiming to be efficient on and portable across exascale high-performance computers. It is motivated by the need to develop a SDM that can model warm-cloud microphysics in domains large enough to resolve shallow mesoscale cloud organisation 𝒪(100 km). This paper introduces Cleo's fundamental structure. First, it explains the choices we made to exploit shared memory spaces on single nodes of high-performance computers. In particular, we discuss our reasoning behind Cleo's ordered and contiguous arrays of structures for superdroplets and grid-boxes, and how we can allocate resources economically and portably through Cleo's MPI domain decomposition, Kokkos library thread parallelism, and coupling to a host dynamical driver. Second, this paper explains how Cleo uses monoids to construct arbitrary combinations and permutations of microphysical processes and data output. These monoids make microphysics and data output highly and easily configurable, which facilitates warm-cloud process understanding. Furthermore, they reduce conditional code branching and enhance Cleo's readability and maintainability. Performance tests on single nodes of Levante HPC System show Cleo has optimal linear scaling with increasing number of superdroplets as well as sufficient memory- and strong-scaling on both CPUs and GPUs to invite CLEO's next stage of development, its optimisation for large eddy simulations on multiple nodes of top-tier and exascale high-performance computers.

1 Introduction

Persistent and large discrepancies between observations of the spatial and temporal patterns of warm-rain and Large Eddy Simulations (LESs) have been ascribed to cloud microphysics for decades (e.g. Randall et al., 2003; Abel and Shipway, 2007; Ackerman et al., 2009; vanZanten et al., 2011) (see also summaries in Khain et al., 2015; Morrison et al., 2020). With the advent of Global Storm Resolving Models (GSRMs) which remove the need to parametrise convection, cloud microphysics is now one of the leading sources of uncertainty in global models too (e.g. Miyakawa et al., 2014; Stevens et al., 2020; Suematsu et al., 2021; Bao and Windmiller, 2021; Lang et al., 2023; Takasuka et al., 2024; Naumann et al., 2025). Much of the uncertainty in both LES and GSRMs arises from epistemic uncertainty, that is uncertainty caused by our insufficient knowledge of microphysical processes and condensate attributes. However, distinct from that, these models also suffer from uncertainty caused by the Eulerian microphysics schemes we conventionally use to represent the knowledge we already have (Morrison et al., 2020). To tackle uncertainty in both LES and GSRMs we not only need to address fundamental gaps in our knowledge, but we must also reconsider the way we represent cloud microphysics in models (Grabowski et al., 2019).

Conventional, Eulerian bulk/moment and bin microphysics schemes have intrinsic deficiencies, for example due to numerical diffusion, the categorisation of condensates into discrete types, and gross assumptions about particle size distributions. The consequent deficiencies in simulated clouds are well documented and contribute to uncertainties in their radiative properties as well as precipitation (e.g. Stevens and Seifert, 2008; Khain et al., 2015; Jian et al., 2021; Schulz and Stevens, 2023). To state just one example, numerical diffusion in physical space is known to activate droplets at cloud-edge and broaden the droplet size distribution, and so may cause unphysical precipitation from stratocumulus (Stevens et al., 1996; Hoffmann, 2016; Morrison et al., 2018; Dziekan et al., 2019). Problems like these are inherent to bulk and bin microphysics schemes and therefore cannot be eradicated by refinement of such models. Indeed, decades of research has aimed to limit the impact of numerical diffusion on cloud microphysics, and yet the problem still persists (Grabowski et al., 2019). If we seek to resolve such issues definitively, then the fundamental limitations of bulk and bin schemes need to be surmounted.

The Super-Droplet Method (SDM; Shima et al., 2009), is a fundamentally different model of cloud microphysics with different limitations and a number of key conceptual advantages (Grabowski et al., 2019). In SDM the cloud condensate population is modelled by “superdroplets”. Superdroplets' motion, microphysics and attributes (for example their masses and radii) are similar to what we expect of real condensate particles; however a superdroplet also has a “multiplicity”, which expresses how many real particles that superdroplet represents. SDM is a much closer representation of the underlying physics than bulk and bin schemes and, because it is a Lagrangian model, it is non-diffusive (in both physical and mass space) (Morrison et al., 2018; Grabowski et al., 2019). Several SDM implementations already exist, (e.g. Shima et al., 2009; Arabas et al., 2015; Naumann and Seifert, 2015; Brdar and Seifert, 2018; Jaruga and Pawlowska, 2018; Shima et al., 2020; Chandrakar et al., 2021; Bartman et al., 2022; de Jong et al., 2023; Matsushima et al., 2023) and are similar to other particle-in-cell-based cloud microphysics models (e.g. Sölch and Kärcher, 2010; Andrejczuk et al., 2010; Riechelmann et al., 2012; Hoffmann et al., 2015), the key difference being SDM's profound convergence property. As the number of superdroplets increases, the model, including its algorithm for droplet collision-coalescence, tends towards what we expect from explicit individual particle simulations (Shima et al., 2009). Although this does not eliminate the epistemic uncertainty caused by insufficient knowledge of microphysical processes and condensate attributes, SDM's convergence property makes its physical interpretation both conceptually elegant and straightforward. This in turn helps us explore such knowledge gaps.

SDM is also particularly well-suited to trends in high performance computing. Aside from converging, its collision-coalescence algorithm conserves the number of simulated particles, i.e. superdroplets, which is beneficial for load balancing. Moreover SDM is embarrassingly parallelisable because its algorithms for microphysics and motion act independently on individual superdroplets or non-overlapping pairs of superdroplets, which makes it ideal for multi-threaded CPU and GPU architectures and explains impressive scaling with increasing number of superdroplets (Arabas et al., 2015; Bartman and Arabas, 2021; Dziekan and Zmijewski, 2022; Matsushima et al., 2023). Its scaling with increasing microphysical complexity is also better than that of bin models once the number of superdroplet attributes is greater than about four (Shima et al., 2009). This makes SDM better suited to modelling droplets with more complexity, which is rarely done in bin and bulk schemes and is necessary to better represent important microphysical processes, for example riming in mixed-phase clouds (Brdar and Seifert, 2018) and aerosol chemistry in aerosol-cloud interactions (Jaruga and Pawlowska, 2018). Both the conceptual and computational advantages of SDM make it a promising tool for better modelling of cloud microphysics (Morrison et al., 2020), and in step with the growth of high performance computers (HPCs) we can now use SDM for more ambitious simulations than ever before.

One emerging application for SDM is to model warm-rain in domains large enough to resolve shallow mesoscale organisation 𝒪(100 km). The evaporation of warm-rain drives cold pool formation and is therefore a key player in mesoscale circulations and cloud organisation (Barnes and Garstang, 1982; Seifert, 2008; Seifert and Heus, 2013; Vogel et al., 2021). Of particular importance in the tropics, mesoscale cloud organisation has major repercussions for the radiative properties of shallow clouds (Nuijens and Siebesma, 2019; Bony et al., 2020) and is highly sensitive to the formation of precipitation (Seifert and Heus, 2013; Yin et al., 2024). There is also ample evidence that mesoscale cloud organisation modulates warm-rain formation (e.g. Stevens et al., 2005; Nuijens et al., 2009; Schulz et al., 2021; Radtke et al., 2022, 2023), and so there are important two-way interactions between the micro- and meso-scale of warm-clouds which, in order to be understood, require an accurate representation of the droplet size distribution, a clear depiction of the microphysical processes at play, and domains large enough to resolve shallow mesoscale circulations. SDM would be an ideal candidate to meet such requirements, except that, presumably limited by computational resources, the largest domains for SDM have not yet exceeded 𝒪(10 km) (Sato et al., 2017; Chandrakar et al., 2021; Matsushima et al., 2023; Yin et al., 2024) – far smaller than the 𝒪(100 km) domains needed for resolving mesoscale organisation, and far smaller than that which may be feasible with the growth of exascale high-performance computers, computing at least 10¹⁸ floating-point operations per second. This is just one example of how, if SDM could be made inexpensive enough at large-scale, it would benefit the community who use large-domain explicitly-resolved-dynamics models – including not just LESs, but also GSRMs.

Hence we are creating Cleo: a novel implementation of SDM aiming to model warm-cloud microphysics efficiently on top-tier and exascale high-performance computers. Cleo's name is a tribute to two “super-women”, Cleopatra VII Thea Philopator and the pseudonymous mathematician Cleo.

In this first paper we introduce Cleo's fundamental computational structure and explain how it is intended to facilitate warm-cloud process understanding and to exploit single nodes of high-performance computers. Specifically, this paper explains how we exploit the massively parallel heterogeneous resources and the hierarchical memory layout of individual nodes of high-performance computers. Notwithstanding SDM's many advantages in comparison with bulk and bin microphysics, SDM faces the typical challenges of Lagrangian particle models; namely simulations can be memory-intensive and it is time-consuming to transport particles in memory. The particle transport's performance on high-performance computers depends on both its shared- and distributed-memory parallelism, and the distributed-memory parallelism strongly varies not only with Cleo's fundamental structure, but also with the fluid-flow given by the host dynamical driver as well as the computer network topology. In a follow-up paper we will therefore present Cleo's distributed-memory parallelism (MPI domain decomposition) and evaluate its performance specifically for LES across multiple nodes of top-tier and exascale high-performance computers. With regard to Cleo's fundamental computational structure on single nodes, the memory and particle transport expenses are prevalent when accessing, sorting, and shuffling superdroplets. Conscious of this, we chose Cleo's memory layout with the aim of economically allocating memory and optimising memory access patterns for loops over superdroplets, whilst compromising on other aspects of performance, namely Single Instruction, Multiple Data (SIMD) parallel processing and scalability with increasing number of superdroplet attributes. Secondly, we have sought to maximise Cleo's freedom to allocate resources on massively parallel heterogeneous computer architectures; both through its handling of domain decomposition, which allows SDM and the dynamics to distribute work independently of one-another, as well as through its use of Kokkos for portable thread parallelism (Trott et al., 2022, 2021), in doing so we aim to be able to allocate resources efficiently.

As well as for high performance, Cleo has various features to facilitate warm-cloud process understanding. We took the idea of monoid sets from mathematics and created a computational representation for them which enables adaptive time-stepping, avoids conditional code branching, and makes microphysics and data output highly configurable. The flexibility of microphysics makes sensitivity studies easy to conduct and, for data output, it's easy to target output specifically to ones needs. Also for the purpose of sensitivity studies we made it straightforward to switch Cleo's grid and couple Cleo to different dynamical drivers. The companion to this paper, Bayley et al. (2026), describes the numerical methods for the warm-cloud microphysics and droplet motion which Cleo makes available.

Cleo's fundamental design is described in Sects. 2–4, with particular emphasis on how it is intended to make SDM simulations suitable for high performance computers (HPCs). In Sect. 2 we describe Cleo's underlying memory layout and how it is designed for efficient memory management; Sect. 3 is for how Cleo enables economical resource allocation, and Sect. 4 describes Cleo's flexibility through monoids. Overviews of Cleo's time-stepping routine and performance on a single node are found in Sects. 5 and 6, respectively.

2 Memory Layout

Efficient memory management is principally achieved by Cleo's underlying memory layout of superdroplets and grid-boxes (Sect. 2.1 and 2.2). We store superdroplets within each shared memory space in an ordered and contiguous chunk of memory. This optimises access patterns for loops over superdroplets and avoids needing additional buffers which increase memory consumption. However, this layout comes at the cost of sorting superdroplets and increasing cache-misses as the number of superdroplet attributes increases. We explain the details behind this compromise in Sect. 2.3.

2.1 Grid-Boxes

“Grid-Boxes” (cells) compose Cleo's spatial domain. Each grid-box defines a distinct region of space for the SDM collision algorithm and thereby determines which superdroplets may interact with one another. In some sense the volume of each grid-box therefore determines the accuracy of the SDM because as the volume of each grid-box decreases given a certain number of superdroplets, the resolution of the flow-field for superdroplet motion increases, and the multiplicity of each superdroplet decreases, thus decreasing the errors in the SDM.

Grid-Boxes are structures consisting of three parts as shown in Fig. 1a. The “GBx-index” is a unique immutable identifier, typically a constant unsigned integer; the “state” defines the macrophysical properties of the volume, such as the wind velocity and thermodynamics (temperature, pressure, relative humidity etc.); and the “view” of superdroplets specifies all the superdroplets which occupy the volume at a given time.

Figure 1 Schematic for the memory layout of (a) each grid-box, and (b) each superdroplet. The grey box contains a key for the memory consumption of each unit.

The state may change due to microphysical processes, for example condensation lowers the water vapour pressure and increases the temperature, or it may change through information received from the host dynamical driver as explained in Sect. 3.3.

The view specifies all the superdroplets at a given time whose spatial coordinates fall within the grid-box's boundaries. Rather than actually containing the superdroplets, it is an indicator to their location(s) elsewhere in memory. For example, the view could be a linked-list of pointers to individual superdroplets, or it could be two pointers to the start and end of a sub-section of a larger array of superdroplets (as it is in our current implementation).

Grid-Boxes are agnostic to the underlying grid, meaning they contain no information about their spatial coordinates or their neighbours. This is because LES grids can vary substantially, for example SCALE uses a Cartesian Arakawa-C grid (Sato et al., 2015; Nishizawa et al., 2015), whereas ICON uses a icosahedral-triangular Arakawa-C grid (Hohenegger et al., 2023), and decoupling the microphysics from the grid of the dynamics driver makes it simple to switch between different grids to aid model inter-comparison studies. For a given grid, a collection of maps/functions called “grid-maps” are assembled, to map from a given GBx-index to grid-dependent information. Such maps take a given GBx-index and return for example the coordinates of the boundaries of that grid-box in a particular direction or the GBx-index of the neighbouring grid-box. In doing so they determine the type of grid, the boundary conditions of the model, and whether the model is 0-D, 1-D, 2-D, or 3-D. Thus, changing these properties for different simulations only requires a replacement of the grid-maps. The rest of Cleo, including the grid-boxes themselves, remains unchanged.

2.2 Superdroplets

Each superdroplet in Cleo is a structure composed of four, or optionally five parts, as shown in Fig. 1b. The four compulsory parts are its attributes, multiplicity, spatial coordinates, and “SD-GBx-index”. The SD-GBx-index always matches the GBx-index of the grid-box whose boundaries enclose the superdroplet's spatial coordinates. It is not a requisite of SDM but assists with keeping the superdroplets ordered by grid-box, which is necessary for efficient cache loading as explained in Sect. 2.3.

Optionally, each superdroplet can be given an unique identity to facilitate tracking the evolution of individual superdroplets. Alternatively, when tracking is not required, superdroplet identities can be given no unique address. During compilation these identities are optimised out, enabling a simple memory-saving measure especially for simulations involving a very large number of superdroplets. The same optimisation could be applied to grid-boxes' GBx-indexes (since their positions in memory could act as unique identifiers instead) however we keep GBx-indexes for generality.

Figure 2 Schematic for the AoS memory layout of grid-boxes and superdroplets. For illustration here the coordinates of the first four superdroplets are within the first grid-box, hence their SD-GBx-indexes match that grid-box's GBx-index and that grid-box's view refers to them.

2.3 Access Patterns

The memory layout of both grid-boxes and superdroplets in a single shared memory space targets efficient memory access patterns. In this layout, the structures for each grid-box are stored together in a single array (an Array of Structures; AoS), and likewise superdroplets are stored in their own AoS, in a contiguous block of memory completely separate from grid-boxes. Additionally the superdroplets' array is ordered by increasing SD-GBx-index, and hence the view of superdroplets for each grid-box is formed by references to the start and end of the sub-block of superdroplets whose SD-GBx-indexes match the grid-box's GBx-index. This is all summarised in Fig. 2. The memory layout we have chosen benefits some computational operations at the expense of others. In particular, ordering the superdroplets by SD-GBx-index makes cache loading efficient but requires the additional expense of sorting them during superdroplet motion, and using AoS favours data transfer over SIMD operations. In the following section we discuss in more detail how these conflicts manifest in the algorithms Cleo uses to conduct SDM microphysics.

Single, contiguous arrays of grid-boxes and superdroplets efficiently handle memory allocations. As opposed to having superdroplets scattered in memory, e.g. by having a separate array of superdroplets for each grid-box, having a single array of superdroplets requires less frequent memory (de)allocations when superdroplets move around the domain, and therefore also avoids the need for buffers which would increase the memory consumption of the model.

Figure 3 Schematic of the superdroplet array before, during and after superdroplet motion. To illustrate the ordering, the superdroplets have been given SD-GBx-indexes for grid-boxes from 0 to 9.

The layout we have chosen is also computationally efficient for cache loading during loops over grid-boxes and superdroplets – exactly as occurs during microphysics and while updating superdroplets' spatial coordinates – because both grid-boxes and superdroplets are not only contiguous in memory but also ordered. One advantage is that we maximise locality of reference, in contrast to unordered or non-contiguous data. This means that for loops over grid-boxes, cache lines are more fully utilised because the block of data loaded in the caches contains more of the superdroplets required for each grid-box. This applies even during superdroplet collisions, since in SDM only superdroplets in the same grid-box may interact with one another. For loops over superdroplets, also when they are nested within loops over grid-boxes, maximising locality of reference means that the same grid-box is accessed repeatedly in as fast succession as possible. The second major advantage of our ordered and contiguous memory layout is that during parallelised loops over grid-boxes, threads on the same core which act on neighbouring grid-boxes are more likely to access superdroplets in the same data blocks, whereas threads on different cores are more likely to access independent superdroplet data blocks. Both this and data locality result in fewer cache misses and less contention in memory access, which is particularly important for high-performance computers because they are characterised by large cache lines and numerous threads and cores working in parallel. In many respects our memory layout amounts to cache blocking, which can help to exploit the embarrassingly parallel nature of SDM, as was shown by Matsushima et al. (2023), who employed cache blocking, albeit differently, to drastically improve the performance of another SDM implementation. As such, Cleo is designed for both efficient memory allocation and access patterns for SDM on high-performance computers.

Whilst keeping the superdroplet array ordered by SD-GBx-index makes loops over superdroplets highly efficient, it also means that during superdroplet motion we must sort the array, as illustrated in Fig. 3. Superdroplets' spatial coordinates and hence SD-GBx-indexes may change, breaking the crucial condition for cache blocking that superdroplets occupying the same grid-box exist in a contiguous chunk of memory. We therefore apply a counting-sort algorithm to the superdroplet array as the last step of superdroplet motion. Although this algorithm scales linearly with increasing number of superdroplets, it nevertheless involves time-consuming data copies, scatter-gather patterns and/or atomic operations, and it doubles the memory consumption of the superdroplet array. Choosing to keep the superdroplet array ordered is therefore a trade-off. On the one hand it makes loops over superdroplets more efficient, on the other hand it makes superdroplet motion more expensive due to sorting.

Using an AoS for superdroplets rather than a Structure of Arrays (SoA) prioritises efficient data movement over vectorisation. With the AoS layout, irrespective of the number of attributes of each superdroplet, the speed of algorithms which sort/shuffle superdroplets during motion/microphysics scale only with the number of superdroplets. In contrast any SoA layout which still maintains the advantages for memory access patterns of ordered superdroplet data would need to perform sorting/shuffling on the array for each individual sub-component of the superdroplets separately, resulting in more computations. Another strength of AoS over SoA is that the addition or removal of superdroplets during motion involves less memory (de)allocations and more efficient cache loading. However, with the AoS layout, operations that act on a single attribute can be less efficient because SIMD parallel processing is harder to apply. Also if the number of attributes of each superdroplet is increased, the number of cache misses likewise increases.

Which data layout results in better performance certainly depends on the computer architecture (e.g. cache size and types of processors) as well as the number of superdroplets and number of superdroplet attributes as a result of the compromise between fast computations through vectorisation and fast memory access patterns. For simulations involving many superdroplets (lots of data transfer) with few attributes (few cache misses), the memory layout we have chosen should be more suitable, but rigorous benchmarking would be needed to confirm the preferred layout for a given setup. In the end, some compromise between the two (an Array of Structures of Arrays, AoSoA) may be the better than either extreme. The decision for the memory layout of Cleo is nevertheless bolstered by the work of Matsushima et al. (2023), who also focused on data movement to optimise their SDM's performance, as well as by the rising weight given to the High Performance Conjugate Gradients Benchmark when ranking HPCs, which demonstrates the increasing concern for the performance of memory bandwidth limited computations.

3 Resource Allocation

Cleo has the potential to allocate resources economically because of its parallelism and the way it couples to a host dynamical driver. We used Kokkos to implement thread parallelism and Message Passing Interface (MPI) for domain decomposition. As such, Cleo can make use of the available memory hierarchies and execution resources on a diverse set of HPC architectures. Whereas Cleo advects superdroplets according to the fluid flow (wind fields) itself, to advect thermodynamic fields (temperature, pressure etc.) Cleo must be coupled to a host dynamical driver capable of advection. However the domain decomposition and computational resources for the dynamics can be independent from those used by SDM. This allows us greater freedom to improve load balancing and therefore to run economical simulations.

3.1 Kokkos Thread Parallelism

Cleo depends on Kokkos for performance portable thread parallelism across a diverse range of computer architectures. Kokkos is a collection of libraries which enable us to write a single source code that builds for an arbitrary number of CPU cores and optionally GPU cores (Edwards et al., 2014). During compilation on a given computer architecture, Kokkos evaluates the hardware available and optional user-defined specifications, to determine how to appropriately allocate the memory and execution resources to achieve high performance. For example, it determines whether to use column- or row-major layouts for arrays of grid-boxes and superdroplets in order to avoid cache misses, and it chooses task sizes for threads of parallelised loops which optimise cache blocking. By doing so, Kokkos implements (hierarchical) parallelism of loops over grid-boxes and superdroplets and selects memory layouts designed to be efficient for a specific processor or accelerator. (Alternatively the memory layouts and task distribution can be manually specified.) When loops over superdroplets are nested inside loops over grid-boxes, as occurs during microphysics and superdroplet motion but not necessarily during initialisation or data output, hierarchical parallelism is invoked; that is both the outer loop over grid boxes and the inner loop over superdroplets are parallelised. The exception is if Kokkos (or manually the user) determines that for the given hardware available the number of grid-boxes or superdroplets is too small for it to be worthwhile to parallelise the loop.

The key benefits of using Kokkos are that in abstracting the parallelism, we are able to separate the software's purpose from the details of a computer's hardware. Thus from one source code we can provide specialised thread parallelism at the same time as portability and maintainability across many architectures. At the time of writing, we have tested CUDA, OpenMP and C++Threads, and the full list of available parallel execution spaces also includes OpenMPTarget, HIP, SYCL and HPX.

3.2 MPI Domain Decomposition

For distributed memory parallelism Cleo uses MPI domain decomposition. The layout described in Sect. 2 is replicated for the sub-domain given to each MPI process and an extra step in superdroplet motion between the second and third steps of Fig. 3 is activated to send superdroplets between processes. This is the only MPI communication Cleo requires for its domain decomposition during time-stepping. To choose the decomposition, Cleo first calculates all the possible ways grid-boxes could be distributed across MPI processes such the number of grid-boxes on each process is as equal as possible, then Cleo chooses the domain decomposition from these options which keeps the most vertical columns on the same MPI process. Since the most frequent motion of superdroplets is usually vertically, by choosing this domain decomposition we typically avoid excessive communication compared to the other options which also evenly distribute the number of grid-boxes in the domain.

This is the first implementation of MPI parallelism in Cleo and is certainly not the most optimal. Since the workload on each MPI process actually scales to first order with the number of superdroplets, weighting the domain decomposition by the number of superdroplets in each grid-box, rather than the number of grid-boxes, would be the first improvement one could make to achieve better load balancing. Further refinement could also determine the domain decomposition by estimating the activity of the superdroplets too, based on the idea that superdroplets predicted to be in or near cloudy regions of the domain would generally be expected to collide more frequently and grow/shrink more rapidly via condensation/evaporation, and hence require more floating-point operations than superdroplets in non-cloudy regions. The major challenge with such superdroplet-based domain decompositions is that superdroplets (and cloudy regions) move around the domain and superdroplets may be added or removed depending on a simulation's boundary conditions (e.g. its aerosol/superdroplet sources and lateral superdroplet sinks). The consequent changes in the number and activity of superdroplets changes the workload of MPI processes and a better optimised domain decomposition would dynamically adapt to restore balance. One well-established way of dynamical load-balancing in Lagrangian particle-based methods is using Hilbert space-filling curves to calculate the most evenly balanced partitioning of the total workload across MPI processes. In Cleo, any superdroplet-based domain decomposition would need to obey the added constraint that superdroplets in the same grid-box of the domain remain contiguously ordered in memory, in order to benefit the memory access patterns on each node, as discussed in Sect. 2.3.

The performance of any optimisation of Cleo's MPI domain decomposition strongly depends on the fluid flow given by the host dynamical driver as well as the computer network topology. The fluid flow determines the degree to which superdroplet movement unbalances the load on each MPI process for a given domain decomposition. The computer network topology determines the speed (cost) of MPI communication, which influences the optimal size of the workload given to each process and the acceptable degree of load imbalance for efficient resource allocation. We therefore reserve optimisation of Cleo's MPI domain decomposition and its performance evaluation for a follow-up paper which discusses them in the context of realistic flows and coupling to a specific dynamical driver on top-tier and exascale high-performance computers.

3.3 Coupling Cleo to a Host Dynamical Driver

Cleo is designed to run concurrently to a host dynamical driver called a “dynamics-solver” and exchange information with it via a “dynamics-coupler”. In a one-way coupling, Cleo receives information for the state of each grid-box from the dynamics-solver, whereas in a two-way coupling, Cleo also sends state information back. For illustration, a simple one-way coupling is when thermodynamics at a particular time-step are received from a dynamics-solver which stores arrays read from a file. Given the information stored in the state of each grid-box (wind velocity and thermodynamics), Cleo can enact microphysical processes on the superdroplets in each grid-box and it can move superdroplets according to Fig. 3: first updating their coordinates and SD-GBX-Indexes and then re-ordering the superdroplet arrays. The numerical methods for microphysics and for updating superdroplet coordinates are provided in Bayley et al. (2026). Whilst Cleo can thus advect superdroplets, to advect thermodynamic fields as per the fluid flow, Cleo must be coupled to a dynamics-solver with a fluid-dynamical core.

The dynamics-coupler determines the degree of independence of the grid and of the resource allocation for Cleo and the dynamics-solver. For the maximum possible flexibility, Cleo can be coupled to a dynamics-solver through Yet Another Coupler (YAC; Hanke et al., 2016). YAC uses MPI communication and can also interpolate variables between different grids meaning that, by using YAC, Cleo and the dynamics-solver can not only have different domain decompositions, but also have different grids. Using different grids is advantageous because it enables the dynamics-solver to compose the domain in the optimal way for its fluid-dynamics. Meanwhile Cleo can compose the domain more favourably for SDM, for example using a nested grid, or grid boundaries which reduce grid-box volumes and/or simplify the numerics of superdroplet motion. In general, using MPI means Cleo and the dynamics-solver do not have shared memory and that their domain decompositions are independent. Whilst this can result in costlier communication, it also maximises our freedom to optimise the load balancing and so economise the allocation of computer resources.

4 Flexibility through Monoids

Cleo uses monoids to enable model flexibility without added run-time from conditional code branching. In abstract mathematics, a monoid is a closed set with an associative binary operation and an identity element. Put simply, monoid sets are composed of elements (members) which result in another element of the set when they're added pairwise together, and the identity element is defined such that the outcome of its addition with another element is just that other element unchanged. Acrylic paint is one example of a monoid set, whereby different colours are the elements of the set and transparent acrylic paint is the identity element. Mixing two acrylic paints together is their associative binary operation (in this special case, also commutative).

We apply the idea of monoids computationally by creating objects in code that can be combined pairwise and associatively with one-another. To do this, first we define a set of constraints a code-object must satisfy in order to be allowed to be combined with another such object. In other words, we define constraints an object must satisfy to be part of a certain code-based monoid set. Secondly, alongside these constraints, we define a function (an associative binary operation) which determines exactly how two objects which obey these constraints are combined together such that the return of the function is a new object which itself satisfies the given constraints. In this way, any number and permutation of such objects can be combined sequentially to create the final object that is used at run-time.

As an example, consider some code which calls a function “do_microphysics” through an object called “Z”. Z can be any object as long as it obeys some constraints: it contains a function called “do_microphysics”, which takes superdroplets as an argument and potentially modifies them. Now let us create any number of objects, A, B, C, … which all obey these constraints but have different implementations of the “do_microphysics” function. For example A.do_microphysics may model condensation, whereas B.do_microphysics may model collision-coalescence. We then define an associative binary operation for these objects that returns a new one which likewise satisfies the constraints of the “do_microphysics” function. For example the operation could return “AB” from A and B where AB.do_microphysics = A.do_microphysics followed by B.do_microphysics. Note that since monoids are not necessarily commutative, “AB” is allowed to cause different microphysical outcomes to those of “BA” – in this example, the growth of superdroplets by condensation before collisions can result in different outcomes to collisions followed by condensation (ı.e AB≠BA). Since we chose to define the outcome of the associative binary operation (AB) to also satisfy the constraints of the “do_microphysics” function, we can then combine AB with further objects, C, D, etc. to create one final object, Z, whose “do_microphysics” function is an assembly of all the microphysics we desire (again, associatively but not necessarily commutativity). Everywhere else in the code we use Z, for example to call the “do_microphysics” function (rather than A, B, or C etc.), meaning that if we wish to change the microphysics we need only recompile and run the model with a different combination of things comprising Z. This allows us extraordinary freedom regarding the implementation of the “do_microphysics” function whilst keeping the rest of the code intact.

To implement monoids in Cleo, we take advantage of template meta-programming with C++20 concepts. Specifically, for a given monoid set we use a C++20 concept to impose the constraints of the set on templated types, and we overload the right-shift operator (operator≫) to implement the function for the associative binary operation. We then create structures/classes which obey the C++ concept and can thus be combined during compilation, in whatever way desired, to instantiate the templated type that is used throughout the code at run-time. Using C++20 concepts is not essential, but has the added benefit of more comprehensible error messages than ordinary template meta-programming in C++.

There are many other plausible ways to express monoids in code. Within older C++ standards (pre-2020) one could instead use for example a base class, a trait, or even just constraints written on a piece of paper. However, base classes would create more restricted, closed definitions of monoids, and the latter options would be more verbose and error-prone. Based on how we defined computational monoids above, other compiled or just-in-time compiled languages could also implement monoids computationally with the same overall behaviour as our monoids in C++. Naturally, languages with object-oriented programming paradigms are better-suited to implementing monoids since monoids are fundamentally described by set theory.

Using monoids to enable model flexibility provides faster run-time performance than using conditional code branching. The decisions regarding which functions are called occurs upon instantiation of each templated type during compilation. This allows compilers to efficiently inline and optimise the used code, whilst at the same time increasing code locality by not compiling unused code. More importantly, the monoidal code is in stark contrast to a code which makes decisions at run-time, for example by using if-statements to determine whether or not to call functions. Indeed, demanding the same degree of flexibility from a code which uses if-statements would require an excessive number of conditional branches – in the example above, to be able to model just A, B, AB, or BA requires three. If this were implemented in Cleo, each branching would add computational expense and reduce the code's modularity. Modular as opposed to monolithic code is more readable and easier to maintain because it can be understood in piecewise units and each unit can be tested and developed independently, as well as concurrently. Although at the time of writing Cleo does not yet have unit tests for its monoids, they are still modular units of code and thus facilitate code readability and maintainability.

Thus, to ensure run-time performance whilst still enabling maintainability and flexibility, Cleo defines monoids for microphysics and data output. Section 4.1 details our monoid for microphysical processes, namely the C++ concept which defines microphysical processes, the structure which acts as the identity element of the set, and the right-shift operator used to combine microphysical processes with one another. Likewise, Sect. 4.2 details the monoid set for “observers”, the elements of which can be used for data output.

Figure 4 (a) The flow diagram for the time-stepping routine of some microphysical process “C” as described by the example in the text. (b) The flow diagram between C.run_step and C.next_step when $C=A\oplus B$.

4.1 Microphysical Processes

Microphysics in Cleo is enacted by the instantiation of a templated type which satisfies the constraints of a “microphysical process”. The constraints are chosen to ensure a type instantiated to perform microphysics will function in the time-stepping routine displayed in Fig. 4a. They are that every microphysical process has three (GPU compatible) functions, next_step, on_step, and run_step, each with specific signatures as imposed by the C++20 concept shown in Listing 1. Naturally, types which count as valid microphysical processes may contain additional information beyond these three functions, e.g. a time-step value or a random number generator for a specific microphysics algorithm, however such additional information is not a requirement of a microphysical process.

Listing 1 The definition of a microphysical process in Cleo.

We also choose to define the associative binary operation $A\oplus B=C$ for types that obey the constraints of a microphysical process such that:

• C.next_step = minimum(A.next_step, B.next_step),

• C.on_step = (A.on_step or B.on_step),

• C.run_step = A.run_step, then B.run_step.

This definition enables us to compile a single microphysical process that obeys the time-stepping routine in Fig. 4a but is actually constructed from the sum of any permutation of microphysical processes, as demonstrated in Listing 2. The flow diagram between C.run_step and C.next_step in Fig. 4a is simply supplanted by calling the run_step function of each of the original microphysical processes sequentially, as exampled by Fig. 4b.

Listing 2 Example code showing how microphysical processes are combined in Cleo. Here we choose “≫” to be the operator in the definition of $A\oplus B=C$ for microphysical processes.

By defining the addition of microphysical processes in this way, we also ensure adaptive time-stepping in which the combined microphysical process always takes the largest possible time-step such that the time-steps of original processes are respected. Consider the microphysical process $C=A\oplus B$. The time-stepping routine in Fig. 4a starts at time = 0 s, which is less than the final time for time-stepping, t_final. When C.run_step is first called, both A.on_step and B.on_step return false and afterwards C.next_step updates the time to the smaller value out of A.next_step and B.next_step. On the next loop iteration, the microphysical process(es) out of A and B which return true from their own on_step function enact microphysics, C.run_step then returns and time is incremented again. The loop for microphysics enacts this logic repeatedly until t_final is exceeded and the simulation terminates. For demonstration, if t_final = 20 s and A and B have constant time-steps of 3 and 5 s, respectively, time advances with the letter in the bracket indicating which microphysics would be enacted at each time-step in the following sequence: 0 s → 3 s (A) → 5 s (B) → 6 s (A) → 9 s (A) → 10 s (B) → 12 s (A) → 15 s (A and B) → 18 s (A) → 20 s (B). In general, A and B need not have constant time-steps and can themselves be formed from the summation of other microphysical processes. Thus we can construct a single microphysical process during compilation which is the sum of any permutation of microphysical processes and enables adaptive time-stepping which adheres to each process' time-steps.

Figure 5 (a) The flow diagram for Cleo's time-stepping routine. Each step is coloured by which of Cleo's structures it involves. The steps involved in the SDM sub-time-stepping routine are shown in (b) and are an extension of the logic presented in Fig. 4a to include output and superdroplet motion as well as microphysics.

For completeness of the monoid for microphysical processes, the set has an identity element, NullMicrophysicalProcess shown in Listing 3. This is useful as an initial microphysics construction, or in simulations when no microphysics is desired.

Listing 3 The identity element of the microphysical process monoid in Cleo.

4.2 Observers

Data output, for example writing information to non-volatile memory or statements to a terminal window, is entirely analogous to microphysics. Output is done by the instantiation of a templated type which satisfies the constraints of an “observer” – a concept for things which have the possibility to read and copy information from Cleo at selected points during time-stepping. Structures called “Monitors” are used by some observers for output which cannot be obtained at the start or end of a time-step, but rather must be monitored during a SDM time-step – for example the mass change due to condensation/evaporation. Monitors are also a monoid in Cleo, but we omit the set's definition here for the sake of brevity. The constraints on an observer are chosen to ensure the type instantiated as the observer will function in Cleo's time-stepping routine shown in Fig. 5 (pink boxes). They are that every observer has 6 functions, before_timestepping, after_timestepping, next_step, on_step, at_start_step, and get_monitor, each with specific signatures as imposed by the C++20 concept shown in Listing 4.

Listing 4 The definition of an observer in Cleo.

Any type which satisfies these constraints is a valid element of the monoid for observers and can therefore be combined with another element using the associative binary operation $A\oplus B=C$, which we choose to define as:

• C.before_timestepping = A.before_timestepping, then B.before_timestepping,

• C.after_timestepping = A.after_timestepping, then B.after_timestepping,

• C.next_step = minimum(A.next_step, B.next_step),

• C.on_step = (A.on_step or B.on_step),

• C.at_start_step = A.at_start_step, then B.at_start_step,

• C.get_monitor = A.get_monitor ⊕ B.get_monitor,

as demonstrated in Listing 5.

Listing 5 Example code showing how observers are combined in Cleo. Here we choose “≫” to be the operator in the definition of $A\oplus B=C$ for observers. Monitors are also monoid (also with their own C++20 concept, associative operation, and identity element) but we omit their definition here for brevity.

Exactly like, but independent of, microphysical processes, the sum of different observers ensures adaptive time-stepping. Likewise this means observers can have different output frequencies and ones that are not necessarily factors of each other. In fact, the constraints on the time-stepping functions for observers and microphysics are so similar they could stem from the same parent time-stepping concept, but for simplicity we keep them distinct.

Observers that write data to non-volatile memory are designed to write large datasets efficiently. Buffers hold output data in RAM until they reach a specified chunk size (∼ 10 MB by default). These are then written to chunks of arrays in a dataset according to the Zarr storage specification version 2. Zarr is chosen because it is well-suited to high throughput of large N-dimensional arrays, especially because it enables parallelised reading/writing. The dataset is compliant with the requirements of Xarray, and is therefore also compatible with NetCDF. Data can be written as any data-type of the NetCDF Climate and Forecast (CF) Metadata Convention, but to save memory it is cast to 4-byte types by default. Additionally, arrays of superdroplet data (e.g. their spatial coordinates) are output according to the CF Metadata Convention's contiguous ragged representation. This is because the number of superdroplets may change over time if superdroplets leave or enter the domain and so the length of data at each output time-step may change. A ragged array, as opposed to a fixed shape array with dimensions of the largest length of data, does not contain empty data points and so saves memory.

For completeness, the monoid for observers has an identity element, NullObserver, shown in Listing 6. This is useful as an initial observer construction, or in simulations which do not require output from Cleo.

Listing 6 The identity element of the observer monoid in Cleo.

5 Time-stepping Cleo

Cleo's entire time-stepping routine is shown in Fig. 5. In summary, one time-step consists of seven key stages:

1. Determine the next time-step, t_next.

2. If on a coupling time-step, receive thermodynamics and wind velocity fields for SDM from the dynamics-solver via the dynamics-coupler.

3. If on an observation time-step, perform observation.

4. Time-step the dynamics-solver from the current time to t_next.

5. Time-step SDM from the current time to t_next via the sub-time-stepping routine shown in Fig. 5b, enacting microphysics followed by superdroplet motion.

6. If on a coupling time-step, send thermodynamics from SDM to the dynamics-solver via the dynamics-coupler.

7. Increment the time to t_next.

The length of one time-step and the exact implementation of each of these parts is determined by types instantiated for the dynamics-solver, dynamics-coupler, grid-maps, observer, superdroplet motion, and microphysical process. For example, when using a one-way coupling the dynamics-coupler's function for Step 6 is empty. The sub-time-stepping routine for SDM shown in Fig. 5b is simply an extension of the logic from Fig. 4a to include observers and superdroplet motion as well as microphysics. Although currently it is performed afterwards, to reduce run-time the entire sub-time-stepping routine could in principle be made asynchronous to the time-step of the dynamics-solver.

Table 1 Configuration and software of Levante HPC for the performance tests of Cleo on a single node with only CPUs (for Serial, OpenMP, and C++Threads builds), and with CUDA enabled.

6 Performance

Here we show the performance of Cleo's superdroplet-scaling and strong-scaling using OpenMP, C++Threads and CUDA parallelism. The simulations are performed on a single node of Levante HPC at the German Climate Computing Center (DKRZ) using 1, 16, 64 or 128 CPU threads and, when CUDA is enabled, an additional GPU¹¹. The node's configuration and the software used are listed in Table 1. To change the problem size we decrease the volume of each grid-box, holding the domain volume and the number of superdroplets per grid-box fixed (Matsushima et al., 2023) (as opposed to fixing the grid-box volumes and increasing the number of superdroplets per grid-box; Dziekan and Zmijewski, 2022). This not only decreases the variance in the number of droplet collisions and increases the precision of the droplet size distribution, but also increases the resolution of the flow field and the accuracy of each superdroplet's position (Bayley et al., 2026). The better this superdroplet-scaling, i.e. grid-box scaling, the more efficiently we can run accurate simulations.

Figure 6 An example of the performance tests with 16 384 grid-boxes. (a–d) The initial conditions, (e) The growth and trajectories of a random sample of 500 superdroplets between 1500 m < x < 3000 m up-to 80 min for the test case with 4.194 × 10⁶ superdroplets overall.

The simulations are of a 3-D Cartesian domain 6000 m × 300 m × 1500 m with 256 superdroplets per grid-box, varying grid-spacing as listed in Table 2 and with the initial conditions shown in Fig. 6a–d. They include condensation, evaporation, collision-coalescence and droplet motion using the configuration listed in Table 3 for the physics described in Bayley et al. (2026). The droplets in every grid-box are initialised from a bimodal log-normal dry-aerosol distribution of NaCl. The modes are at 20 and 150 nm with standard deviations of 1.4 and 1.6, respectively, their relative probability is 3:2, and the total droplet number concentration is 500 cm⁻³. Superdroplets are created by randomly sampling 256 evenly-in-log-space bins between 5 nm and 1 µm – as in the “Single-SIP-init” method of Unterstrasser et al. (2017), but in radius not mass space. Superdroplets' spatial coordinates are uniformly randomly sampled from the space within each grid-box. The dynamics are a simplified version of the 2-D kinematic flow model from Arabas et al. (2015). To extend the simulation into 3-D, we replicate the 2-D (x-z) flow uniformly along the additional (y) dimension. We prescribe the same divergence-free wind field and hydrostatic equilibrium, but the kinetics are held constant in time, without feedback from microphysics nor relaxation towards the initial profiles. To mimic sub-cloud and in-cloud conditions, the temperature has a lapse rate of 9.8 K km⁻¹ below 750 m and 6.5 K km⁻¹ above, while the water vapour has a lapse rate of 2.97 $\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ below 750 m and super-saturation is fixed at 0.1 % above. The surface pressure is 1013.15 hPa, the surface temperature is 297.9 K, and the surface water vapour mass mixing ratio is 16 g kg⁻¹. We ran the simulations for 80 min without data output and present results per 1 s of simulated time. An example of the superdroplets' evolution up-to the end of the simulation is shown in Fig. 6e.

Table 2 The problem size for the performance tests in terms of the number of grid-boxes composing the domain and the consequent total number of superdroplets and domain grid-spacing.

Rogers et al. (1993)Simmel et al. (2002)

Table 3 SDM configuration for the tests of Cleo's superdroplet- and strong-scaling.

Figure 7 The superdroplet-scaling of the wall-clock time for each simulated second given the maximum set of resources for each build on Levante HPC: (a) Serial (dash-dotted), (b) CUDA (dotted), (c) OpenMP (solid), and (d) C++Threads (dashed). The colours decompose the total time (blue) into the time spent on SDM (cyan), motion (purple) and microphysics (red).

The superdroplet-scaling of the wall-clock time is shown in Fig. 7 for the maximal set of resources for each build: serial (1 thread), CUDA (128 CPU threads and 1 GPU), OpenMP (128 threads) and C++Threads (128 threads). The wall-clock time for all the builds is almost entirely time-stepping SDM, as opposed to initialisation or the dynamics-solver, and time-stepping SDM scales linearly with increasing number of superdroplets once the total number of superdroplets is large enough. This is because the bottleneck of microphysics is the random shuffling of superdroplets required by the method for SDM collisions, and the bottleneck of motion is the aforementioned counting-sort algorithm (Sect. 2.3). Both of these algorithms scale linearly with increasing number of superdroplets. For C++Threads and OpenMP builds, motion (i.e. sorting) is the dominant bottleneck, while for serial and CUDA builds, it is the microphysics (i.e. shuffling). The time-per-call of motion is, however, three times greater than shown here since the time-step was 3 s. The expense of motion per-call is therefore comparable to that of microphysics for all builds, and considerably more for OpenMP and C++Threads.

Figure 8 The strong-scaling of the speed-up of wall-clock time relative to the serial simulations with increasing number of CPU threads on Levante HPC. The line-style denotes the build and the colours from purple to green indicate increasing problem size in terms of the total number of superdroplets in the domain. The speed-up of motion from CUDA exceeds the y-axis scale and is reported in Table 4 instead.

The strong-scaling in Fig. 8 shows the speed-up of the wall-clock time with increasing number of CPU threads for various problem sizes¹. Overall, none of the builds reach optimal efficiency (Fig. 8a); CUDA speed-up is limited by microphysics (Fig. 8b), while OpenMP and C++Threads speed-ups are limited by motion (Fig. 8c).

The limit for CUDA parallelism due to microphysics can be seen in Table 4. The speed-up of motion maximizes at 150 for the largest problem sizes, however microphysics appears to plateau at smaller problem sizes and at 70. The bottleneck is again the shuffling algorithm we use during collision-coalescence, presumably due to the relatively slow clock-speed of a single GPU thread during the serial step of the Fisher-Yates algorithm we use for shuffling. A different (possibly parallelised) algorithm may therefore improve performance, or alternatively, the shuffling could be performed concurrently with other parts of the time-stepping routine. For OpenMP, the strong-scaling of microphysics is almost perfect, at 100 % efficiency for large enough problem sizes, however, the overall speed-up is limited by the speed-up of motion which is no greater than 8. The scaling of C++Threads is better than OpenMP despite slightly less-optimal scaling of the microphysics because motion scales better for all but the largest problem size. The bottleneck of motion on CPUs is not the data copying itself, but rather the synchronisation of threads required by scatter-gather patterns in the counting-sort algorithm. Data copying nevertheless doubles the memory footprint of the superdroplets. A more specialised sorting algorithm which exchanges rather than copies superdroplets and uses the relationship between the GBx-indexes to avoid scatter-gather patterns would therefore provide better performance both in terms of memory usage and speed. Our results are in agreement with Matsushima et al. (2023) who found that after optimising their SDM, superdroplet tracking (movement and sorting) was indeed the bottleneck of SDM on CPUs, and for the largest simulation this was predominately due to superdroplet sorting.

Table 4 The speed-up of wall-clock time for the CUDA build relative to the serial simulations for problem sizes plotted in Fig. 8.

The bottlenecks in Cleo's microphysics and motion are expected to become less acute with more complex simulations. The addition of more microphysical processes, for example to model frozen condensates, will increase the cost of microphysics in ways that decrease importance of the shuffling algorithm. Coupling Cleo to a dynamics-solver capable of advection will diminish SDM's role in the total run-time and introduce other dynamics-related bottlenecks. Above all, the expense of transporting superdroplets between MPI processes is expected to make the cost of sorting during superdroplet motion inferior to the cost of MPI domain decomposition, even once optimised.

Figure 9 (a) The superdroplet-scaling of the serial build for the resident set size high-water memory consumption (brown) and for the maximum memory allocation (grey); (b, c) The strong-scaling of the same memory footprints for the parallelised builds: OpenMP (solid), C++Threads (dashed) and CUDA (dotted). The colours from purple to green indicate increasing problem size in terms of the total number of superdroplets in the domain.

Table 5 The apparent size of superdroplets and grid-boxes calculated from their contribution to the high-water memory consumption for the simulations with 32 768 grid-boxes (256 superdroplets per grid-box).

Figure 9 shows the superdroplet- and strong-scaling of Cleo's memory consumption is as expected. The memory footprints increase linearly with increasing number of superdroplets and are almost constant with increasing CPU threads. The overhead from using OpenMP or C++Threads parallelism is almost insignificant. It is larger for smaller problem sizes but still no greater than 1 %, and is predominately from the creation of thread-safe parallel random number generators and arrays for the scatter-gather patterns. Due to the large number of threads on a GPU, the overhead of CUDA parallelism for problem sizes less than approximately 10⁶ superdroplets can be relatively large because thread-safe parallel random number generators are more expensive than the superdroplets and grid-boxes themselves. For large enough problem sizes however, the overhead of CUDA is also insignificant, and the maximum memory allocation (i.e. the superdroplet memory allocation) can actually be substantially smaller than in serial. The high-water memory consumption, that is the largest resident set size memory used, can be divided into the apparent memory consumed by a single superdroplet and grid-box. As demonstrated in Table 5, for all builds both superdroplets and grid-boxes are effectively 𝒪(100 bytes), whereas by default stand-alone they are 68 and 120 bytes, respectively.

7 Conclusions

To satisfy our curiosity about cloud responses to climate change, SDM offers a promising new way to model cloud microphysics. SDM simulations, even in the absence of collision-coalescence, overcome many of the problematic features of conventional microphysics representations and this has already proven fruitful, for example in studying the role of turbulence during droplet growth (Li et al., 2018, 2020; Grabowski and Thomas, 2021; Chandrakar et al., 2024; La et al., 2025). However, the computational cost currently inhibits SDM from certain areas of research, an important case being the study of interactions between cloud microphysics and their mesoscale organisation. This research requires LES that can clearly depict warm-cloud microphysical processes and resolve droplet sizes whilst also capturing mesoscale circulations. We therefore need a SDM for studying warm-clouds in LES with hectometre resolution and 𝒪(100 km) domains. Cleo is a new computational implementation of SDM intending to be fit for that purpose. Of course, to make such large-domain LES feasible, Cleo must be able to harness both the shared- and distributed-memory parallelism of exascale computers, and whilst the shared-memory performance depends mostly on the fundamental structure of the code, the distributed-memory parallelism is also sensitive to the fluid-flow given by the host dynamical driver and the computer network topology. As such, we will evaluate Cleo's distributed-memory parallelism in an follow-up paper specific to LES conducted across multiple nodes of top-tier and exascale high-performance computers. This paper meanwhile discusses Cleo's fundamental computational design and it's performance on shared-memory architectures (single nodes). In Bayley et al. (2026), the numerics of Cleo's warm-cloud microphysics is described.

Efficient memory management is a principal feature of Cleo's design for high computational performance of large simulations. We have chosen the layout of grid-boxes and superdroplets in contiguous ordered chunks of memory to try to minimise cache loading and contention, and support data locality. Whilst keeping superdroplets ordered requires a costly sorting algorithm, this expense will become less significant with larger simulations involving MPI domain decomposition, more complex microphysics and a fluid-dynamical core. The memory layout we have chosen also favours data transfer over SIMD operations, as is deemed appropriate for simulations involving frequent transport of a large number of superdroplets which have a low number of attributes. Since simulations of many superdroplets are exceptionally memory demanding, Cleo's optional features, such as superdroplet tracking, can be optimised out of simulations. Data output is also designed economically by following the Zarr storage specification to write binary files and by using ragged arrays to represent superdroplet data.

Cleo's resource allocation is designed flexibly so that large simulations can be economised. To exploit parallelism on HPCs, we use MPI for domain decomposition and Kokkos for abstracted and hierarchical thread parallelism. Kokkos makes Cleo concise and maintainable, but most importantly performance portable, meaning Cleo can make use of the various forms of thread-parallelism available on a diverse set of computer architectures, including those with GPUs. Further flexibility is built into how Cleo couples to a host dynamical driver. The coupling allows SDM and the dynamics to have independent domain decompositions, meaning we have greater freedom to allocate computational resources in a way that achieves load balancing and therefore produces more economical simulations.

Cleo is highly configurable as well as computationally performant. Monoids, which we define through template meta-programming with C++20 concepts, enable adaptive time-stepping and make choices regarding which variables to output and which microphysical processes to include highly flexible. In contrast with flexibility obtained through conditional code branching, flexibility via monoids has faster run-time and makes the code more readable. The coupling between Cleo and a host dynamical driver not only allows different domain decompositions, but also makes it easy to switch host model to conduct sensitivity studies and optimise the numerical methods of SDM and the dynamics-solver separately. Additionally, Cleo has flexibility over the choice of grid for SDM. All this flexibility makes Cleo well-suited for detailed study of microphysics, for example evaporation, its influence on hydrometeor evolution and atmospheric dynamics.

Cleo's current state of development raises opportunity for further optimisation and larger simulations. The bottlenecks of Cleo's algorithms for microphysics on a single node are the shuffling during SDM collisions (acute in GPU simulations), and the sorting during superdroplet motion (in CPU-only simulations). Using a different shuffling algorithm may be more optimal for GPUs, but optimisation could also come from making parts of the shuffling algorithm execute concurrently with other computations at each time-step. Similar optimisations could provide better performance of the sorting algorithm on CPUs, as well as lower its memory consumption. However, whether such optimisations are worthwhile should be assessed with more ambitious simulations, where the transport of superdroplets between nodes, additional microphysics, and/or the expense of the dynamics-solver and communication with it, could surpass these algorithms as the major bottlenecks. Given Cleo's current state of development however, we can speculate about the minimum cost of such ambitious simulations. Consider a LES with a 150 km × 150 km × 7 km domain with grid-spacing 150 m horizontally, 40 m vertically and 128 superdroplets per grid-box. In other words, 𝒪(10¹⁰) superdroplets overall. For a conservative estimate on the cost, let us assume that the overhead from particle transport via MPI and communication with the dynamics-solver increases Cleo's wall-clock by a factor of 10. Then, by using 1000 NVIDIA A100 GPUs (250 “Levante-like” GPU-nodes), the SDM setup from the performance tests presented here suggest that 𝒪(10⁷) superdroplets per MPI process would need approximately one second of wall-clock time for every second of simulated time and have an order of magnitude less memory consumption than the maximum on each node. That would make Cleo's throughput approximately a factor of 10 less than the strong-scaling limit of ICON (Klocke et al., 2025), and even for an extremely conservative estimate, where Cleo's performance is slower by a factor of 100, a simulation of several hours that takes several weeks to run is still conceivable. Future development of Cleo is therefore a promising avenue to pursue towards using SDM to study the interactions between cloud microphysics and mesoscale cloud organisation.