| 30 Jul 2026
GPU-accelerated atmospheric chemical kinetics using single precision in the ECHAM/MESSy (EMAC) model v2.55 with MEDINA v2.0
Abstract. Atmospheric chemical kinetics puts a major computational burden on Earth system models, yet its embarrassingly parallel nature of the problem structure makes it well suited to GPU acceleration. Here we present a single-precision floating point (FP32) GPU implementation of the MECCA chemical kinetics submodule in the ECHAM/MESSy Atmospheric Chemistry (EMAC) model. It is enabled by the MEDINA code-to-code compiler v2.0 and compatible with the Rosenbrock-family numerical solvers generated by the Kinetic Pre-Processor (KPP). The approach exploits FP32 arithmetic inside the GPU kernels while preserving double precision in the host EMAC workflow via explicit type conversion on the GPU.
We evaluate the numerical fidelity and achieved performance in both an idealised box-model configuration and in full EMAC production simulations on the JUWELS Booster system. Across tested mechanisms and resolutions, FP32 reproduces the double-precision CPU reference to within <0.004 % in total-species-mass differences, with discrepancies confined to chemically depleted tracers near machine representation limits. In production-scale runs, FP32 further reduces the average chemical-kinetics execution time by about 30 % relative to the FP64 GPU implementation and lowers the memory footprint of FP-intensive data structures. Chemistry calculations in FP32 enhance scalability for high-resolution and chemically complex simulations.
The manuscript presents a useful, in-scope engineering and validation study extending MEDINA to single-precision GPU chemistry, with solid accuracy (<0.004% mass error) and performance (~30-45% kernel speedup) results across both box-model and full production EMAC configurations. I recommend the manuscript for publication subject to minor revisions; addressing the suggestions below would further strengthen the paper.
Comments
Consider adding how the chemistry-kernel speedup translates into whole-model wall-clock savings
The comparison between CPU-side and GPU-side typecasting (Sec. 2.1) is stated only qualitatively ("extensively testing... found to require overall more execution time"). Please consider adding quantitative numbers.
Line 22: Please consider explaining why adaptive time-step control is relevant enough to justify the phrase "the numerical integration of the stiff chemical ODEs remains computationally demanding despite recent optimization of the adaptive time-step control."
Since floating-point error can accumulate nonlinearly over longer integrations, consider discussing whether the reported <0.004% mass differences are expected to remain stable over longer runs, or over time during the one-month simulation.
Line 135: Consider adding the comparison of the memory footprint against approach 1, for example: requires ∼1.5× the memory footprint of a double-precision-only GPU kernel on the device (∼3× that of Approach 1's device allocation, since Approach 1 only stores the single-precision array on-device).
Figure 2 compares FP32-GPU against FP64-CPU, which conflates precision differences with CPU/GPU architectural differences. The authors do isolate architectural noise via an FP64-CPU-vs-FP64-GPU control run (L178-179, ≤0.0004%) and do provide a same-architecture FP32-GPU vs. FP64-GPU comparison for the full EMAC model (Fig. 5), but not for the box-model case shown in Fig. 2. For consistency with Sec. 3.2 and to fully isolate the precision effect in the box-model results as well, consider adding a box-model panel analogous to Fig. 5 (FP32-GPU vs. FP64-GPU), or at minimum reporting the FP64-CPU-vs-FP64-GPU control run visually alongside Fig. 2.
Table 2 reports speedup as a percentage using S[%] = (t_DP − t_SP)/t_SP × 100 (Eq. 1), normalized by the new (SP) runtime rather than the baseline (DP) runtime — a valid but less standard convention that risks being misread as (t_DP − t_SP)/t_DP. For consistency with the "Speedup ×" notation used in Fig. 8 (Eq. 4), and for more intuitive interpretation, consider reporting Table 2 as a dimensionless ratio (t_DP/t_SP) instead of a percentage. Note this would not make the two speedup metrics identical — Table 2 compares SP vs. DP at fixed solver/resolution, while Fig. 8's Si is normalized to a fixed single-node DP baseline for scaling — but using the same "×" format throughout would avoid the current mix of two different speedup conventions (Eq. 1 vs. Eq. 4) in the same paper.
Table 2 shows a speedup of only ~42–45% from FP64 to FP32, well below the ~2× that might naively be expected either from halving data volume (memory-bound case) or from the GPU's ~2× FP32:FP64 FLOP throughput ratio (Table 1: 19.5 vs. 9.7 TFLOPS). Table 3 shows the memory-traffic reduction itself is non-uniform across the cache hierarchy (20.75% at L1-L2, 63.24% at L2-DRAM), suggesting the kernel's speedup is not simply explained by a uniform halving of memory traffic. Since the Introduction motivates FP32 partly via the Roofline model (Williams et al., 2009), consider adding a brief discussion — ideally informed by a roofline-style analysis or profiling breakdown — of whether the Rosenbrock integration kernel is memory- or compute-bound, and whether non-floating-point operations (integer indexing, pivoting, adaptive step-size control/branching) that don't benefit from reduced precision are capping the achievable speedup below the theoretical 2×.
Line 235 / Figure 7: It is unclear what generates the "distribution" of execution times shown for each node/GPU configuration in Figure 7. For instance, does each box/violin represent the set of individual per-timestep integration times (ts) accumulated over the course of a single 31-day production run (consistent with Eq. 3, AET = mean of ts over Nsteps)? The text introduces Eq. 3 immediately after this sentence but never explicitly links it to what is plotted in Fig. 7. Please state explicitly, either in the text or the Figure 7 caption, the source and sample size so that the spread shown in the figure can be correctly interpreted.
Consider adding a test that fixes the hardware (number of nodes/GPUs) while varying the problem size — e.g., sweeping the box-model vector length parameter (VL_GLO), which the paper already uses for this purpose in Sec. 3.1 (fixed at 576,000). This would show whether the FP32 speedup reported in Table 2 is already saturated at the tested workload size or would continue to grow with more data per GPU, which is directly relevant given that the paper cites Kerkweg et al. (2025) finding higher speedups (5.14×–7.69×) for larger vectors of chemical species (L96-97). Establishing this relationship for the FP32-vs-FP64 comparison specifically would strengthen the scalability claims made in Sec. 3.2.1 and Fig. 8.