Enhancing the advection module performance in the EPICC-Model V1.6.0 via GPU-HADVPPM4HIP V1.0 coupling and GPU-optimized strategies

The rapid development of Graphics Processing Units (GPUs) has established new computational paradigms for enhancing air quality modeling efficiency. In this study, the heterogeneous-compute interface for portability (HIP) was implemented to parallel computing of the piecewise parabolic method (PPM) advection solver (HADVPPM) on China's domestic GPU-like accelerators (GPU-like), resulting in a GPU-accelerated version denoted as GPU-HADVPPM4HIP V1.0. Computational performance was enhanced through three strategic optimizations: reducing the central processing unit (CPU) and GPU (CPU-GPU) data transfer frequency, thread-block coordinated indexing, and the Message Passing Interface (MPI) and HIP (“MPI+HIP”) hybrid parallelization across heterogeneous computing clusters. Following validation of the GPU-HADVPPM4HIP V1.0 program's offline computational consistency and the pollutant simulation performance of the Emission and atmospheric Processes Integrated and Coupled Community version 1.6.0 (EPICC-Model V1.6.0) on the Earth System Numerical Simulation Facility (EarthLab), comprehensive performance testing was conducted. Offline benchmark results demonstrated that, for a one-dimensional input array containing 10⁸ double-precision floating-point values, GPU-HADVPPM4HIP V1.0 compiled with the −O3 optimization option achieved a maximum speedup of 556.5× on a GPU-like accelerator relative to the Fortran HADVPPM implementation compiled with the −O0 option as the baseline. Integrating GPU-HADVPPM4HIP V1.0 into EPICC-Model V1.6.0 yielded three distinct versions: the initial HIP-based version (HIP-Ori), a version optimized for CPU and GPU communication frequency (HIP-Opt1), and a further-optimized version employing a thread-block coordinated indexing strategy (HIP-Opt2). Compared to the HIP-Ori version, HIP-Opt1 achieved a model-level computational efficiency improvement of 17.0×. Building upon HIP-Opt1, HIP-Opt2 delivered an additional 1.5× enhancement in computational efficiency. At the module level, including CPU and GPU data transfer overhead, the GPU implementation improves computational efficiency of the advection module by 39.3 %; when communication cost is excluded, the advection module attains a 20.5× acceleration relative to its CPU counterpart. This coupling establishes a foundational framework for adapting air quality models to GPU-like architectures and identifies critical optimization pathways. Moreover, the methodology provides essential technical support for achieving full-model GPU implementation of the EPICC-Model, addressing both current computational constraints and future demands for high-resolution air quality simulations.

1 Introduction

Air pollution, including fine particulate matter in both urban and rural areas, is associated with an elevated risk of strokes, heart diseases, lung cancer, and acute and chronic respiratory diseases (Atkinson et al., 2010; Kim et al., 2015; Liu et al., 2016; Milton and White, 2020). The air quality forecasting system centered on the air quality model plays a critical role in the timely dissemination of forecasting alerts and early warning information to the public. The accuracy of air quality forecasting is jointly constrained by the spatial resolution of input datasets, including emission inventories, terrain, and meteorological parameters (Gupta and Mohan, 2015; Georgiou et al., 2022). High-resolution model configurations have demonstrated improvements in the accuracy of air quality forecasting (Georgiou et al., 2018; Podrascanin, 2019; Adani et al., 2022; Gao and Zhou, 2024). However, current operational forecasting systems predominantly employ horizontal resolutions ranging from several to tens of kilometers (Wu et al., 2014; Guevara et al., 2021; Tang et al., 2022; Gao and Zhou, 2024), which inadequately address the requirements for urban-scale high-resolution forecasting and precision management.

Computational demands emerge as a critical limiting factor for high-resolution air quality modeling. On the one hand, doubling the horizontal resolution quadruples the number of computational grids. On the other hand, maintaining numerical integration stability necessitates proportional reduction in temporal integration steps (Georgiou et al., 2022). These combined effects result in exponential growth of computational workload with increasing resolution. It is estimated that when the horizontal resolution of the air quality model is increased by 18 times, the computational load of the model increases by 300 times (Thompson and Selin, 2012).

Enhancement of computational efficiency in air quality modeling has been predominantly achieved through hardware-based acceleration strategies. Wang et al. (2017) ported the Global Nested Air Quality Prediction Modeling System (GNAQPMS) to the second-generation Intel Xeon Phi processor, Knights Landing (KNL), achieving a 3.5× computational speedup through hybrid Message Passing Interface (MPI) and Open Multi-Processing (OpenMP) parallelization, vectorization, memory-access optimization, reduction of thread-local storage, and global communication optimization. The gas-phase chemistry module is widely recognized as the dominant computational bottleneck in air quality models, typically accounting for over 40 % of total simulation time (Elbern, 1997; Linford et al., 2011; Wang et al., 2017; Cao et al., 2023). To address this limitation, Wang et al. (2019) developed the multi-point Carbon Bond Mechanism Z (MP CBM-Z) by implementing vectorized computation techniques within the CBM-Z framework. Leveraging Single Instruction Multiple Data (SIMD) architecture, their approach enabled multi-point parallel computation for gas-phase chemistry, achieving a 4.9× acceleration in the chemistry module and a 2.22× overall speedup for the entire NAQPMS model when deployed on Intel Xeon Gold 6132 CPUs.

In recent years, GPUs have emerged as transformative accelerators in artificial intelligence and high-performance computing, driven by their massive parallel computing capabilities. In December 2024, the 64th TOP500 list ranked the El Capitan system first, with a performance of 1.742 EFLOP s⁻¹ in the High Performance Linpack (HPL) benchmark, which measures the sustained double-precision floating-point performance achieved when solving a dense system of linear equations (TOP500, 2024). This computational supremacy primarily originates from its AMD Instinct MI300A GPU accelerators, each containing 14 592 stream processors and delivering a double-precision floating-point performance of 61.3 TFLOP s⁻¹. Remarkably, the computational efficiency of a single MI300A GPU exceeds 1.8 times the peak performance of the Earth-Simulator supercomputer (CPU-based architecture) in Japan, which was the fastest supercomputer in 2003.

The formidable computational capacity of GPUs has opened new directions for enhancing the computational efficiency of air quality models. Alvanos and Christoudias (2017) developed a software package for the ECHAM/MESSy Atmospheric Chemistry (EMAC) model, enabling the automated generation of Compute Unified Device Architecture (CUDA) kernels for the numerical integration of atmospheric chemical kinetics using the Kinetic PreProcessor (KPP; Damian et al., 2002). Subsequent memory optimization and thread management strategies achieved a 20.4× acceleration for the chemistry module on NVIDIA P100 GPUs. In parallel efforts, Sun et al. (2018) implemented CUDA-based optimization for the second-order Rosenbrock chemical solver (Sandu et al., 1997) within the CAM4-Chem global chemistry-climate model. Through strategic enhancements in fully interleaved memory layout, CUDA streams, and constant memory, they achieved an 11.7× speedup for computation alone and a 3.8× speedup when the data transfer between the CPU and GPU is considered on the NVIDIA Tesla K20X GPU. Notably, Quevedo et al. (2025) adapted the third-order Rosenbrock solver in the Community Multiscale Air Quality (CMAQ) model by converting the Fortran code to CUDA Fortran, evaluating its performance across three chemical mechanisms: RACM2, CB6R5, and SAPRC07. Comparative analysis revealed 51 %, 50 %, and 35 % computational efficiency gains on NVIDIA RTX 2080 Ti GPUs, respectively, while maintaining numerical consistency with CPU-based benchmarks. Through code refactoring from Fortran to standard C and the use of the heterogeneous-compute interface for portability (HIP) programming framework, Cao et al. (2025) successfully parallelized the fourth-order Rosenbrock solver on a GPU-like architecture. Concurrently, the total model elapsed time was reduced by 46.9 %. Regarding the advection module, another computational hotspot in air quality models, Cao et al. (2023, 2024) implemented GPU-accelerated adaptations of the advection module in the Comprehensive Air Quality Model with Extensions (CAMx) using CUDA and HIP, respectively, and the optimized advection module achieved maximum speedups of 80.2× on NVIDIA Tesla V100 GPUs and 28.9× on GPU-like.

Following Cao et al.'s (2025) successful implementation of parallel computing for the gas-phase chemistry module in the EPICC-Model on GPU-like accelerators, the computational time proportion of this module was significantly reduced. Consequently, the advection module has emerged as a computational hotspot with comparable time consumption to the optimized chemistry module. To address this shift, this study focuses on enhancing the advection module performance in EPICC-Model V1.6.0 via GPU-HADVPPM4HIP V1.0 coupling and GPU-optimized strategies. Section 2 details the EPICC-Model's computational framework, baseline performance tests, and the heterogeneous computing platform employed in this research. Section 3 elaborates on the optimization framework specifically designed for the computational characteristics of the EPICC-Model advection module. Section 4 presents experimental results, including simulation performance and computational performance analysis.

2 The EPICC-Model and experiments

2.1 The framework of the EPICC-Model

The emission and atmospheric processes integrated and coupled community model version V1.6.0 (EPICC-Model V1.6.0; EPICC-Model Working Group, 2026a; Wang et al., 2025) is an air quality modeling system specifically designed for complex air pollution in China (Zhu et al., 2023), and developed by the Institute of Atmospheric Physics, Chinese Academy of Sciences based on the Earth System Numerical Simulation Facility (EarthLab, Chai et al., 2021). The model framework is fundamentally based on the following species continuity equation and is used to simulate the complex physical and chemical processes of pollutants in the atmosphere.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\partial C_i}{\partial t}}+u{\displaystyle \frac{\partial C_i}{\partial x}}+v{\displaystyle \frac{\partial C_i}{\partial y}}+\mathit{\omega }{\displaystyle \frac{\partial C_i}{\partial z}}\\ & ={\displaystyle \frac{\partial }{\partial x}}\left(K_x{\displaystyle \frac{\partial C_i}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(K_y{\displaystyle \frac{\partial C_i}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(K_z{\displaystyle \frac{\partial C_i}{\partial z}}\right)\\ & +S+P-R_{\mathrm{d}}-W_{\mathrm{ash}}\end{array}\end{array}$$

where C_i is the concentration of pollutant species i, u, v, and ω are the wind-velocity components in the x, y, and z directions, respectively. K_x, K_y, and K_z are the turbulent diffusion coefficients in the longitudinal, latitudinal, and vertical directions, respectively. S represents the emission, P represents chemical transformation, R_d represents dry deposition, and W_ash represents wet scavenging.

The EPICC-Model V1.6.0 adopts a modular architecture developed using the Fortran programming language, a high-performance computing language specifically designed for scientific applications. The model code is open-source and shared (EPICC-Model Working Group, 2026b). This open-source code repository enables rapid integration of novel mechanisms and modules proposed by diverse research groups, thereby enhancing collaborative development efficiency.

The computational framework and workflow of the EPICC-Model V1.6.0 are illustrated in Fig. 1. The system primarily comprises three components: model inputs, physical-chemical processes, and outputs analysis. Model input data include emissions, meteorological data, and other datasets. Emission inventories such as the Multi-resolution Emission Inventory for China (MEIC; Li et al., 2017), Emissions Database for Global Atmospheric Research (EDGAR; Crippa et al., 2024), Hemispheric Transport of Air Pollution version 3 (HTAP_v3; Crippa et al., 2023), and Inversed Emission Inventory for Chinese Air Quality (CAQIEI; Kong et al., 2024) can be used. Meteorological data are predominantly derived from simulations generated by the mesoscale Weather Research and Forecasting (WRF) model. Other input datasets include configuration files, terrain data, photolysis data generated by the Tropospheric Ultraviolet and Visible (TUV) radiation model, initial conditions (ICs), and boundary conditions (BCs). Physical-chemical processes primarily include horizontal advection, vertical diffusion, dry deposition, wet scavenging, gas-phase chemistry, aqueous-phase chemistry, heterogeneous reactions, inorganic aerosol thermodynamics, etc. For vertical diffusion, either the scheme of Byun and Dennis (1995) or the Yonsei University (YSU) scheme (Hong et al., 2006) can be selected to calculate the turbulent vertical diffusion coefficient. The dry deposition module can employ either the scheme of Wesely (1989) or Zhang et al. (2003) to compute deposition velocities. The gas-phase chemistry module offers the option to utilize either the CBM-Z (Zaveri and Peters, 1999) or CB6r5 (Yarwood et al., 2020) chemical mechanisms. For heterogeneous reactions, the model defaults to the scheme of Li et al. (2012). Additionally, it integrates mechanisms for HONO heterogeneous chemical reactions (Zhang et al., 2022), sulfate heterogeneous chemical reactions (Li et al., 2018), and N₂O₅ heterogeneous hydrolysis (Yang et al., 2024). Inorganic aerosol is simulated using the ISORROPIA aerosol thermodynamic model (Nenes et al., 1998). The aqueous-phase chemistry module originates from the Regional Acid Deposition Model (RADM, Chang et al., 1987). Regarding model output analysis, the EPICC-Model can generate pollutant concentration fields, pollutant deposition fluxes, process analysis data, and source apportionment data. For a comprehensive technical description of the model architecture and implementation details, refer to the EPICC-Model Working Group (2026a).

Figure 1 Computational framework and workflow of EPICC-Model V1.6.0. In the physical and chemical processes component, yellow denotes physical modules and orange denotes chemical modules. ICs and BCs denote initial conditions and boundary conditions, respectively. FNL denotes the Final Operational Global Analysis produced by the National Centers for Environmental Prediction, GFS denotes the Global Forecast System, and WPS denotes the WRF Preprocessing System, where WRF represents the Weather Research and Forecasting model.

For the horizontal advection module that is the focus in this study, two high-precision numerical schemes are available, the positive-definite mass-conservative differencing scheme (Walcek and Aleksic, 1998) and the piecewise parabolic method (PPM, Colella and Woodward, 1984). The PPM scheme, an extension of high-order Godunov's method, operates by partitioning the integration domain into subregions and approximating solutions using parabolic functions. Renowned for its numerical precision and robustness in complex fluid dynamics, this classic algorithm has been widely adopted in atmospheric chemistry models including the latest CMAQ and CAMx (Appel et al., 2021; Emery et al., 2024). Within the EPICC-Model V1.6.0 framework, the advection module sequentially executes transport processes in the x-direction and y-direction. It employs species-specific PPM solvers (HADVPPM subroutine) for gaseous species, inorganic aerosols, organic aerosols, dust, and sea salt. Our previous studies have demonstrated a significant acceleration performance of PPM solver in the CAMx model through HIP heterogeneous programming technologies for GPU-like (Cao et al., 2024).

2.2 Benchmark performance of testing

As mentioned above, Cao et al. (2025) implemented the fourth-order Rosenbrock solver for the gas-phase chemistry module in the EPICC-Model, employing the CBM-Z chemical mechanism. Through code restructure from Fortran to standard C programming language and implementation of the HIP heterogeneous programming framework, the computational efficiency of the gas-phase chemistry module improved by 2.88 times when accounting for data transfer between CPUs and the GPU-like accelerators.

Figure 2 illustrates the changes in the computational time proportion among modules before and after the heterogeneous parallel implementation of the gas-phase chemistry module on GPU-like accelerators, as reported by Cao et al. (2025). Specifically, Fig. 2a shows the time proportion of each module in the original Fortran version, while Fig. 2b presents the corresponding proportion after the gas-phase chemistry module was ported for parallel computing on GPU-like. As shown in Fig. 2, the implementation of parallel computing for gas-phase chemistry modules on GPU-like achieved significant efficiency improvements, reducing its computational time proportion from 45.7 % to 29.9 %. Notably, the computational time proportion of MPI_Barrier synchronization function has decreased from 23.9 % to 13.4 %. A critical observation emerged regarding the advection module, whose computational time proportion increased from 13.3 % to 25.0 %, establishing it as a new performance bottleneck comparable to the optimized chemistry module. This performance shift necessitates subsequent optimization efforts focusing on heterogeneous porting and parallel acceleration of the PPM scheme for GPU-like architectures within the EPICC-Model framework, aiming to enhance the computational efficiency of the advection module.

Figure 2 The computational time proportion among modules (a) for the Fortran version and (b) after implementing the gas-phase chemistry module on GPU-like accelerators for parallel computation.

2.3 Hardware platform and software environment of experiments

All performance testing of the EPICC-Model V1.6.0 and the heterogeneous adaptation and optimization of its advection module were conducted at the Earth System Numerical Simulation Facility (EarthLab; Chai et al., 2021). Jointly developed by the Institute of Atmospheric Physics, Chinese Academy of Sciences, and collaborating institutions, EarthLab is designed for Earth system modeling and high-resolution regional environmental simulations and employs a CPU-GPU heterogeneous architecture. Within this platform, the term “GPU-like accelerator” refers to a domestically developed general-purpose accelerator in China that performs the same massively parallel acceleration role as a conventional GPU and supports GPU-style programming through the HIP framework. The GPU-like accelerator used in this study has 3840 computing units and 16 GB of device memory. As summarized in Table 1, each GPU-like computing node contains two of China's domestic CPU processors and two GPU-like accelerators connected through Peripheral Component Interconnect Express (PCIe) 4.0 buses, and both the CPUs and GPU-like accelerators are first-generation devices. Cao et al. (2024) previously used HIP to accelerate the PPM advection solver of the CAMx air quality model on the same type of domestic GPU-like accelerator, while Wang et al. (2021) demonstrated the large-scale application of the HIP-based LICOM3 ocean model. The software environment includes the Intel OneAPI 2021.3.0 toolkit for compiling the CPU code and the dtk-23.04.1 toolkit for compiling the GPU-like code.

Table 1 The hardware components and software environment for the dedicated accelerator node on the EarthLab.

Compared to CPU processors, GPU-like accelerators demonstrate superior capability in launching massive thread-level parallelism. Similar to the NVIDIA GPU architectures (NVIDIA, 2020), these GPU-like accelerators employ a three-level parallelism hierarchy comprising grids, blocks, and threads, which collaboratively execute parallel computations through coordinated indexing. Specifically, a computational grid is partitioned into multiple thread blocks with three-dimensional coordinates, each thread block containing an array of three-dimensional indexed threads. As the fundamental execution unit, individual threads perform concrete computational tasks, each possessing a unique index ID that precisely determines its spatial position within the thread block hierarchy. Consequently, the design of hierarchical indexing schemes coordinating blocks and threads constitutes a critical challenge in achieving efficient parallel computation for three-dimensional numerical modeling grids.

Analogous to AMD's Radeon Open Compute (ROCm) software stack (AMD, 2023), the dtk-23.04.1 toolkit (Cao et al., 2024) includes programming models, tools, compilers, libraries, and runtimes for artificial intelligence and high-performance computing applications on GPU-like accelerators. Mirroring ROCm's design paradigm, dtk-23.04.1 adopts the HIP programming language as its application programming interface (API). This implementation leverages the Single-Instruction Multiple-Thread (SIMT) execution model to effectively manage and coordinate massive thread parallelism on GPU-like accelerators.

3 Implementation details

3.1 Description of the heterogeneous porting and optimization scheme

The heterogeneous porting and parallel optimization schemes of this study are illustrated in Fig. 3. Similar to the heterogeneous porting approach for the advection module of the air quality model CAMx on GPU-like accelerators (Cao et al. 2024), the first step involved the porting and adaptation of the HADVPPM advection solver from the EPICC-Model V1.6.0 to GPU-like accelerators. Firstly, the Fortran code of the HADVPPM subroutine was reconstructed using standard C programming language, followed by implementing parallel computing on GPU-like accelerators through the HIP API. Similar to CUDA program execution on NVIDIA GPUs, the implementation of GPU-HADVPPM4HIP V1.0 on GPU-like accelerators follows four key steps: (1) Device memory allocation via the hipMalloc interface, (2) Data transfer from CPU to GPU-like accelerator through hipMemcpy operations, (3) Parallel computation using kernel launching (hipLaunchKernelGGL) with thread-index-based parallel processing after successful data transmission, and (4) Final data retrieval from GPU back to CPU through hipMemcpy operations.

Figure 3 Heterogeneous porting and parallel optimization scheme of advection module in the EPICC-Model V1.6.0.

Following the implementation of GPU-HADVPPM4HIP V1.0 parallel computing on GPU-like accelerators, the second phase involves architecture-specific parallel optimizations tailored for GPU-like characteristics. Three optimization strategies were sequentially implemented to fully exploit the SIMT vectorization parallelism of GPU-like accelerators, thereby enhancing the computational performance of the EPICC-Model advection module. These strategies include: (1) reducing the frequency of communication between the CPU and GPU, (2) collaborative indexing between threads and blocks, and (3) hybrid parallelization of “MPI+HIP”. For systematic reference, three progressively optimized configurations were designated, namely HIP-Ori, HIP-Opt1, and HIP-Opt2. The HIP-Ori is baseline implementation after GPU-HADVPPM4HIP V1.0 integration into EPICC-Model without optimizations. The HIP-Opt1 is the version implementing reduced CPU-GPU communication frequency. The HIP-Opt2 is the enhanced version incorporating collaborative thread-block indexing. The hybrid “MPI+HIP” parallelization strategy was implemented across all three heterogeneous versions to enhance parallel scalability of the EPICC-Model V1.6.0 on the EarthLab.

3.2 HIP-Opt1: Communication frequency reduction

Influenced by the evolutionary trajectory of high-performance computing, most geoscientific numerical models, including the EPICC-Model, are predominantly coded in Fortran and designed for general-purpose CPU architectures. These models typically execute computations through grid-wise loop iterations. Taking the HADVPPM subroutine in the EPICC-Model as an example, its computational kernel is structured with triple nested loops progressing from innermost to outermost: species loop (loop_species), latitudinal grid loop (loop_j), and vertical grid loop (loop_k). The EPICC-Model innovatively categorizes atmospheric pollutants within the species loop into five distinct classes: gaseous species (74 chemical constituents), inorganic aerosols (14 species subdivided into coarse- and fine-mode particle sizes), organic aerosols, dust aerosols (8 species across 4 size bins), and sea salt aerosols (11 species across 4 size bins).

The EPICC-Model calculates the advection process through looping of chemical-specified variables, which has relatively low computational efficiency. A benchmark test using a two nested grid configuration (horizontal grids: d01 = 228 × 165, d02 = 465 × 300) indicated that each timestep requires approximately 9.8 million calls of the HADVPPM subroutine for advection processes in both x-direction and y-direction. Consequently, the HIP-Ori version, generated by integrating GPU-HADVPPM4HIP V1.0 into the EPICC-Model, incurs 9.8 million CPU-GPU data transfers per timestep. Benchmark tests demonstrated that the computational time for 1 h integration increased from 1015.0 s in the original Fortran version to 20 400.3 s under HIP-Ori version, with frequent CPU-GPU communication identified as a primary performance bottleneck. To address this critical bottleneck, our optimization framework prioritizes architectural redesign of the advection module's loop hierarchy, strategically reducing communication frequency while increasing data transfer sizes to better exploit GPU computational capacity.

Figure 4 illustrates the code-level implementation of CPU-GPU communication optimization, with panel (a) depicting the HIP-Ori baseline and panel (b) presenting the optimized HIP-Opt1 version. In the two nested-domain case, the HIP-Ori configuration required approximately 4.9 million GPU calls for x-direction advection alone. To mitigate this computational overhead, we restructured the advection module's loop architecture and expanded array dimensionality. The HIP-Opt1 optimization framework implements multidimensional array restructuring, beginning with the dimensional expansion of concentration variables from their original 1D representations (Q01d) in HIP-Ori to 4D/5D configurations (Q04d/Q05d), while auxiliary parameters such as grid area adjustment vector (area) and interfacial area adjustment vector (areav) are similarly upgraded from 1D to 3D structures. Prior to GPU execution, these variables undergo systematic multidimensional reorganization, as demonstrated by the transformation of gaseous concentration variables from Q01d_gas(i) to Q04d_gas (i, j, k species) in Fig. 4a–b. This architectural redesign enables complete x-direction advection computation through a single GPU call per pollutant category. Consequently, the total GPU calls for both x and y direction advection decrease from approximately 9.8 million in HIP-Ori to 10 in HIP-Opt1, achieved through one GPU call per spatial dimension across five pollutant categories, thereby optimizing computational efficiency through batched multidimensional data processing.

Figure 4 Code-level implementation of the CPU-GPU communication optimization. Panel (a) shows the baseline Fortran code for HIP-Ori, and (b) shows the optimized Fortran code for HIP-Opt1.

3.3 HIP-Opt2: Thread and block coordinated indexing

The architectural advantage of GPU-like accelerators manifests in their capacity to support massive thread concurrency for parallel computing. To leverage this capability, we implemented the coordinated thread-block indexing method proposed by Cao et al. (2023), in which each grid cell in the two-dimensional horizontal plane is assigned to a dedicated thread. Specifically, blocks were configured based on the meridional grid dimension, with each block allocated threads corresponding to the zonal grid count. This hierarchical parallelization strategy enables parallel processing across the entire two-dimensional horizontal grid through coordinated thread-block resource allocation. To ensure consistency with the original Fortran implementation on the projected grid, the HIP-Opt2 advection routine incorporates the map-scale factors associated with the Lambert projection when calculating horizontal transport. Furthermore, given that GPUs are well-suited for large-scale matrix parallel computations without data dependencies, prior to implementing parallel computation on a two-dimensional grid using thread and block coordinated indexing, we decoupled the iterative computations present in the advection module by introducing intermediate variables. This ensures computational independence at each step, meaning that the computation of the next step does not depend on the results of the previous step, thereby fully leveraging the multi-thread parallel computing capability of the GPU.

3.4 “MPI+HIP” hybrid parallelization

The recent advancements in GPU technology have solidified the dominance of “CPU+GPU” heterogeneous architectures in global high-performance computing, with 9 out of the top 10 supercomputers in 2024 utilizing heterogeneous configurations (TOP500, 2024). These large-scale heterogeneous clusters typically deploy one or multiple GPU accelerators per compute node. Aligned with this trend, EarthLab employs heterogeneous architecture in its dedicated computing nodes, integrating China's domestic CPUs with GPU-like accelerators. To maximize GPU utilization and enhance the parallel scalability of the EPICC-Model on heterogeneous clusters, an “MPI+HIP” hybrid parallelization scheme was designed tailored for EarthLab, inspired by the “MPI+CUDA” approach proposed by Cao et al. (2023) for the CAMx model.

In the original Fortran version of EPICC-Model, the entire simulation domain is decomposed into subdomains using MPI, with each CPU process responsible for computations, including those of the advection module, within one subdomain. In this study, after MPI domain decomposition, we offloaded the advection module computations originally performed on the CPU to China's domestic GPUs through heterogeneous porting, thereby enabling parallel execution of the advection module on GPUs. Considering future high-resolution applications with large-scale data, assigning a single GPU to multiple CPU processes could lead to data-transfer contention and potential GPU memory overflow. Therefore, we adopted an “MPI+HIP” hybrid parallelization approach, in which each participating CPU process is assigned a dedicated GPU accelerator. This design expands the parallel computing capacity of the model on heterogeneous clusters and makes full use of GPU resources. The implementation can be summarized in three main steps: (1) Obtain the MPI process rank information, as well as the number and indices of GPUs available on each compute node; (2) Based on the MPI process rank and using the remainder function in standard C, determine the number and indices of GPUs to be launched; (3) Map each MPI process to a specific GPU index, thereby realizing the “one CPU process – one GPU” configuration in the MPI+HIP hybrid parallel scheme.

As illustrated in Fig. 5, taking the d02 domain configured in Sect. 4.1 with 8 CPU processes and 8 GPU-like accelerators as an example: the EPICC model decomposes the simulation domain into 8 subdomains using the MPI software standard, with each CPU process handling its assigned subdomain. During the execution of the advection module, the “MPI+HIP” hybrid parallel scheme allocates one GPU accelerator to each CPU process. Computational tasks originally performed by the CPU are offloaded to the corresponding GPU; after the advection computations are completed, the results are returned to the CPU.

Figure 5 An example schematic of domain subdivision and mapping to CPU processes, where each CPU process is equipped with a GPU-like accelerator.

4 Experimental results

4.1 Experimental setup

The centre of the simulation domain is located at (35° N, 105° E) and its two true latitude lines are 25 and 45° N, respectively. The EPICC-Model V1.6.0 employs a two-nested configuration, the first domain (d01) covers East Asia with a 45 km × 45 km horizontal resolution on 228×165 grid cells. The lower left corner of the second domain has its starting positions in the grid of the first domain as 36 and 39, respectively, and the second domain focuses on China with a 15 km × 15 km horizontal resolution on 465 × 300 grid cells. In vertical, 20 terrain-following layers are configured with the height of the top layer set to 20 km and six layers below 1 km. The numerical simulation was conducted from 00:00 UTC on 1 July to 23:00 UTC on 31 July 2021, spanning a total duration of 744 h. The initial 168 h period was allocated for model spin-up. The MEIC emission inventory (Li et al., 2017) was adopted as the emission input, and baseline year is 2019. The numerical schemes selected during the EPICC-Model V1.6.0 simulations are listed in Table 2. Furthermore, the PM_2.5 and O₃observations are from the China National Environmental Monitoring Centre, which provides hourly PM_2.5 and O₃ observations for eight cities in China. The station information, including station name and its latitude and longitude, is listed in Table 3. Meteorological inputs are generated using the Weather Research and Forecasting (WRF, Skamarock et al., 2008) model, a mesoscale numerical weather prediction system, and is widely adopted in both theoretical research and operational forecasting. This study employed the WRF version 3.9.1, and the model domain configuration maintains identical nesting architecture and spatial coverage as the EPICC-Model.

Table 2 The physical and chemical numerical schemes selected during EPICC-Model V1.6.0 simulation.

Table 3 The names and latitude-longitude information of the PM_2.5 and O₃ observation stations.

4.2 Simulation performance analysis

4.2.1 Error analysis of the GPU-HADVPPM4HIP

As elaborated in Sect. 3.1, the implementation of the HADVPPM Fortran code on GPU-like comprises two principal phases. Initially, the Fortran code undergoes reconstruction using standard C programming language through a Fortran-to-C conversion process. Subsequently, the standard C-version HADVPPM program is adapted to GPU-like accelerators through C-to-HIP expansion employing the HIP heterogeneous programming technology. Following successful GPU adaptation, the offline precision verification was conducted to compare computational accuracy among three implementations, the original Fortran source code, restructured standard C code, and HIP-accelerated code. To ensure input consistency, a dedicated Fortran program was developed to generate identical input datasets, including 100 double-precision floating-point numbers, for all three implementations. Each implementation executed a complete advection integration computation, with subsequent output recording and analysis. Therefore, the absolute errors (AE) and relative errors (RE) presented in Table 4 represent the average values computed from the 100 double-precision floating-point results produced by the original Fortran code, restructured standard C code, and HIP-accelerated code of the HADVPPM program after performing the advection solution computation under the given input conditions.

Table 4 Comparison of mean AE and RE across compilation configurations for F-to-C, C-to-HIP, and F-to-HIP processes.

Notably, both Fortran and C implementations were executed on China's domestic CPUs and compiled using the Intel OneAPI 2021.3.0 toolkit, while GPU-HADVPPM4HIP was compiled with the dtk-23.04.1 toolkit for GPU execution. For enhanced analytical rigor, we further implemented compilation with −O0 and −O3 optimization flags. The −O0 flag maintains default compilation settings without code optimization, whereas −O3 flag enables more aggressive loop and memory-access optimizations, such as scalar replacement, loop unrolling (Intel Software, 2018).

Table 4 presents the mean absolute error (AE) and relative errors (RE) in computational precision between the Fortran and standard C implementations (F-to-C), standard C and HIP implementations (C-to-HIP), and Fortran and HIP implementations (F-to-HIP) of the HADVPPM program under two compilation configurations. Notably, the −O3 optimization flag prioritizes computational performance through code optimization at the expense of precision degradation. Consequently, the AE and RE values for all three porting processes (F-to-C, C-to-HIP, and F-to-HIP) under the −O0 configuration flag are systematically lower than those under −O3 flag. For instance, between the Fortran and HIP implementations, the AE and RE under −O0 flag are measured at $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{7}}$ % and $\mathrm{4.3}\times {\mathrm{10}}^{-\mathrm{7}}$ %, respectively. However, when compiled with −O3 flag, these errors increase significantly to $\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{7}}$ % and $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{6}}$ %, respectively. Similar error escalation patterns are observed in the F-to-C and C-to-HIP comparisons under −O3 optimization flag. This phenomenon aligns with expectations, as aggressive compiler optimizations (e.g., loop unrolling and memory access reorganization) may introduce numerical instability through altered operation sequences and reduced intermediate precision preservation.

Furthermore, it is noteworthy that across both −O0 and −O3 compilation configurations, the AE and RE values of the F-to-C process consistently exceed those of the C-to-HIP process. This indicates that the computational errors introduced during the heterogeneous porting of the HADVPPM Fortran code from domestic CPUs to GPU-like accelerators predominantly originate from the Fortran-to-C transcoding phase. For instance, under the −O0 configuration, the AE and RE for F-to-C are measured at $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{7}}$ % and $\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{7}}$ %, respectively, whereas those for C-to-HIP are significantly smaller, at $-\mathrm{9.5}\times {\mathrm{10}}^{-\mathrm{9}}$ % and $-\mathrm{8.0}\times {\mathrm{10}}^{-\mathrm{8}}$ %, respectively. This discrepancy arises from inherent differences between Fortran and C in programming paradigms and data precision management. The Fortran-to-C code restructuring involves two fundamentally distinct programming languages, differing in object-oriented design philosophies and numerical representation conventions, thereby introducing computational inaccuracies. In contrast, the HIP programming model, analogous to NVIDIA CUDA, is inherently an extension of standard C. As detailed in Sect. 3.1, HIP achieves GPU compatibility by augmenting standard C programming language with critical GPU-specific functionalities, such as memory allocation and data transfer operations. Since HIP code essentially constitutes an enhanced standard C framework, the C-to-HIP adaptation introduces minimal computational bias, resulting in markedly lower error compared to the F-to-C transformation.

By coupling GPU-HADVPPM4HIP V1.0 with the EPICC-Model and implementing optimizations such as reducing CPU-GPU communication frequency, adopting coordinated thread-block indexing, and employing MPI+HIP hybrid parallelization, we developed the HIP-Opt2 version. Using identical input data and model configurations, the differences between the simulation results produced by HIP-Opt2 and the original Fortran version were evaluated in a real-world case study. The model input data and configurations were set as described in Sect. 4.1. Both versions were integrated for 28 h beginning at 00:00 UTC on 1 July 2021, consistent with the simulation duration used in the performance evaluation in Sect. 4.3.2. The resulting output fields were used for the numerical comparison presented in this section.

Figures 6 and 7 compare the gas and aerosol concentrations calculated by the original Fortran and HIP-Opt2 versions at the end of the 28 h integration. Visually, the two versions produce nearly identical spatial distributions and concentration magnitudes. To reveal the regional distributions of the small numerical differences, species-specific color-bar ranges appropriate for the corresponding difference magnitudes were used in the right-hand panels. The displayed differences are of the order of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$ ppbV for HONO, $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ ppbV for SO₂ and NO₂, and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{2}}$ ppbV for O₃. For the aerosol species, the displayed differences are of the order of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ µg m⁻³ for black carbon (BC) and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ µg m⁻³ for PM_2.5, ASO₄, and ANH₄. Cao et al. (2024) reported that the differences between the HIP and Fortran versions of CAMx were within $\pm \mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{2}}$ ppbV for O₃ and NO₂ and within $\pm \mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{2}}$ µg m⁻³ for ASO₄. By comparison, in the present study, the differences in O₃ and NO₂ are within $\pm \mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{2}}$ ppbV and $\pm \mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ ppbV, respectively, while the difference in ASO₄ is within $\pm \mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ µg m⁻³, all of which are smaller than the corresponding differences reported by Cao et al. (2024). For each species, the grid-cell differences were further averaged. The resulting mean differences are $-\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{9}}$, $\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{7}}$, $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{6}}$, and $\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{4}}$ ppbV for HONO, SO₂, NO₂, and O₃, respectively, and $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{6}}$, $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{5}}$, $\mathrm{4.0}\times {\mathrm{10}}^{-\mathrm{6}}$, and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{6}}$ µg m⁻³ for BC, PM_2.5, ASO₄, and ANH₄, respectively. Taken together, the smaller difference ranges and near zero domain mean differences confirm the close numerical agreement between the Fortran and HIP-Opt2 versions.

Figure 6 Spatial distributions of HONO, SO₂, NO₂, and O₃ concentrations calculated by the original Fortran and HIP-Opt2 versions at the end of the 28 h integration. Panels (a), (d), (g), and (j) show the Fortran results; (b), (e), (h), and (k) show the HIP-Opt2 results; and (c), (f), (i), and (l) show the differences between the two versions. The units of the concentrations and differences are ppbV.

Figure 7 Spatial distributions of BC, PM_2.5, ASO₄, and ANH₄ concentrations calculated by the original Fortran and HIP-Opt2 versions at the end of the 28 h integration. Panels (a), (d), (g), and (j) show the Fortran results; (b), (e), (h), and (k) show the HIP-Opt2 results; and (c), (f), (i), and (l) show the differences between the two versions. The units of the concentrations and differences are µg m⁻³.

All the statistical metrics presented in Table 5 were calculated from the concentration fields at the end of the 28 h Fortran and HIP-Opt2 simulations used to generate Figs. 6 and 7. To further evaluate the suitability of HIP-Opt2 for scientific applications, we followed Wang et al. (2021) and Cao et al. (2024) and calculated the root mean square error (RMSE), standard deviation (SD), RMSE $/$ SD ratio, and normalized mean bias (NMB). RMSE quantifies the differences between HIP-Opt2 and the original Fortran version, while SD represents the spatial variability of each species in the Fortran results. The RMSE $/$ SD ratio therefore measures the magnitude of the computational differences relative to the spatial variability of the corresponding species. NMB measures the overall signed difference of the HIP-Opt2 results relative to the Fortran results. As shown in Table 5, the RMSE $/$ SD ratios range from $\mathrm{6.9}\times {\mathrm{10}}^{-\mathrm{4}}$ % for NO₂ to $\mathrm{9.2}\times {\mathrm{10}}^{-\mathrm{2}}$ % for ANH₄. The NMB values range from $-\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ % for HONO to $\mathrm{2.9}\times {\mathrm{10}}^{-\mathrm{3}}$ % for ASO₄. The small RMSE $/$ SD ratios and near-zero NMB values demonstrate that the computational differences introduced by the HIP implementation are negligible relative to the spatial variability and concentration magnitudes of the evaluated species.

Table 5 Root mean square error (RMSE) between the Fortran and HIP-Opt2 versions, standard deviation (SD) of the Fortran results, RMSE $/$ SD ratio, and normalized mean bias (NMB) of HIP-Opt2 relative to the Fortran version.

4.2.2 Simulation performance verification of the EPICC-Model

After parallelizing the HADVPPM program in the air quality model CAMx on GPU-like accelerators, Cao et al. (2024) conducted comparative analyses of simulation results between NVIDIA GPUs and domestic GPUs through offline and coupled testing approaches. Although both experimental results indicated smaller computational errors introduced by China's domestic GPUs, the study did not validate the discrepancies between CAMx simulation results and actual observational data, particularly regarding model performance in real-case scenarios. To address this gap, the current study integrates GPU-HADVPPM4HIP V1.0 into the EPICC-Model V1.6.0 and performs one-month real-case simulations following the experimental configuration described in Sect. 4.1. This serves dual purposes: firstly, to verify the computational stability of EPICC-Model V1.6.0 in cross-architecture heterogeneous cluster environments using China's domestic hardware, and secondly, to evaluate the model's pollutant simulation performance through observational validation. For observational data, we collected PM_2.5 and O₃ observations from major Chinese cities including Beijing, Taiyuan, Hangzhou, Hefei, Fuzhou, Qingdao, Guangzhou, and Kunming, implementing quality control procedures following the methodology established by Wu et al. (2018). Regarding simulation data, we extracted model outputs from grid cells in the d02 domain corresponding to the geographical coordinates of monitoring stations.

Figures 8 and 9, together with Table 6, present time-series comparisons and statistical evaluations of the simulated and observed daily PM_2.5 and maximum daily 8 h average (MDA8) O₃ concentrations. The formulas used to calculate the statistical metrics are provided in the Supplement. Positive and negative NMB values indicate model overestimation and underestimation, respectively. For daily PM_2.5, the EPICC-Model reproduced the temporal variations reasonably well in most cities. Fuzhou had an RMSE of 3.9 µg m⁻³, an r of 0.87, and an NMB of −14.6 %. Across the eight cities, the mean RMSE, r, and city-level NMB were 9.4 µg m⁻³, 0.70, and 19.7 %, respectively. The positive NMB values in Beijing and Hefei, 89.1 % and 74.0 %, indicate overestimation, whereas the negative NMB values in Guangzhou and Fuzhou, −28.1 % and −14.6 %, indicate underestimation. For MDA8 O₃, the mean RMSE, r, and city-level NMB were 12.1 ppbV, 0.76, and −17.7 %, respectively, indicating an overall tendency toward underestimation. The largest negative NMB values occurred in Taiyuan and Hangzhou, at −41.5 % and −34.3 %, respectively, while Qingdao showed a small positive bias of 4.7 %. Qingdao and Guangzhou had r values of 0.91 and 0.87 and RMSE values of 12.2 and 7.8 ppbV, respectively. These statistical results indicate that the EPICC-Model V1.6.0, integrated with GPU-HADVPPM4HIP V1.0, reasonably reproduces the spatiotemporal characteristics of PM_2.5 and O₃ concentrations on China's domestic heterogeneous computing clusters.

Figure 8 Time series of observed and simulated daily PM_2.5 concentrations in major Chinese cities from 8 to 31 July 2021. Panels (a)–(h) represent Beijing, Taiyuan, Hangzhou, Hefei, Fuzhou, Qingdao, Guangzhou, and Kunming, respectively.

Figure 9 Time series of observed and simulated MDA8 O₃ concentrations in major cities of China from 8 to 31 July 2021. Panels (a)–(h) represent Beijing, Taiyuan, Hangzhou, Hefei, Fuzhou, Qingdao, Guangzhou, and Kunming, respectively.

Table 6 Statistical evaluation of the PM_2.5 and O₃ simulations of the EPICC-Model in eight Chinese cities. NMB denotes normalized mean bias. The average row represents the arithmetic mean of the city-level metrics.

The underestimations or overestimations of PM_2.5 and O₃ simulations observed in specific cities and periods are primarily attributable to two factors. First, the coarse model resolution – the horizontal resolution of the d02 domain in this experiment is 15 km – hindered the accurate representation of terrain features, meteorological variables, and spatial variations in emission sources. Second, discrepancies exist between the baseline year of the emission inventory and the simulation year. Specifically, the MEIC emission inventory used in this study is based on 2019 data, whereas the simulation year is 2021. For cities with stringent pollution control policies, such as Beijing, the 2019 MEIC inventory may inadequately reflect actual emission reductions achieved by 2021, particularly for pollutants under strict abatement measures. This discrepancy could lead to overestimated simulated concentrations.

4.3 Computational performance analysis

4.3.1 Offline performance comparison

A comparison of offline computational results indicates that the computational errors introduced during the Fortran-to-HIP heterogeneous porting process under both compilation settings are minimal. Specifically, the HADVPPM program exhibits small discrepancies on the order of 10⁻⁷ when ported from CPU to GPU-like architectures. That is, differences in the computational results only begin to appear at the seventh decimal place. Based on the verified consistency of offline results, the computational performance of the HADVPPM program was further evaluated on domestic CPU and GPU-like under different compilation configurations. To achieve this, Fortran-based test programs were implemented to generate randomized input arrays with varying scales for both Fortran and HIP versions of the HADVPPM program. These random arrays are one-dimensional and contain 10², 10³, 10⁴, 10⁵, 10⁶, 10⁷, and 10⁸ random numbers. Figure 10 illustrates the computational time and speedup ratios of the HADVPPM program across different compilation options and data scales on domestic CPU and GPU-like. It is explicitly stated that the execution time measurements for the HADVPPM program on the GPU-like exclusively account for the computational time of the kernel functions on the device, while overheads such as GPU memory allocation and host-device data transfer are excluded.

Figure 10 The offline computational time (a) and speedup ratios (b) of the HADVPPM program across different compilation options and data scales on domestic CPU and GPU-like.

As illustrated in Fig. 10a, under both −O0 and −O3 compilation flags, the SIMT vectorized parallel computing advantages of the GPU-like accelerators become prominent for large data scales exceeding 10⁴, demonstrating significantly higher computational efficiency compared to domestic CPU. At a data scale of 10⁸ with the −O0 flag, the Fortran-based HADVPPM program required approximately 11.97 s to complete computation on the CPU, while the HIP version on the GPU-like achieved the same task in 0.55 s, yielding a speedup ratio of 21.87×. The −O3 compilation flag further enhances computational efficiency through automated code optimization, albeit at a slight cost to numerical precision. At the same 10⁸ data scale, the Fortran version on the CPU required 2.75 s, whereas the HIP version on the GPU completed computations in 0.02 s, achieving a remarkable speedup of 128.03×. Notably, the HIP version compiled with −O3 exhibited 556.5× higher efficiency than the Fortran version compiled with −O0 flag. However, for smaller data scales such as 10² or 10³, the architectural advantages of the GPU-like diminish, with computational efficiency comparable to or even lower than that of a CPU. For instance, at a 10² scale with the −O3 option, the GPU's performance matched the CPU (speedup ratio is approximately 1). Under the −O0 option, the GPU's speedup dropped to 0.16×, indicating inferior efficiency relative to the CPU. These results underscore that the GPU-like is well-suited for large-scale matrix parallel computing tasks.

As described in Sect. 3.1, a complete heterogeneous computation of GPU-HADVPPM4HIP V1.0 on the GPU-like involves critical processes such as GPU memory allocation, data transfer between CPU and GPU, and kernel launching for parallel computing. To quantify the overheads of these processes, the time consumption and proportional contributions were analyzed under the −O0 compilation option for different data sizes. In Fig. 11, the labels of “Memory Allocation”, “Memory Copy”, “Kernel”, and “Synchronize” represent the processes of memory allocation on the GPU, bidirectional data transfer between CPU and GPU, kernel launch and parallel computation, as well as thread synchronization within the kernel, respectively.

Figure 11 Time consumption (a) and proportional contributions (b) of GPU memory allocation, CPU-GPU memory copy, kernel execution, and synchronization under the −O0 compilation option for different input-data sizes.

As shown in Fig. 11, for data scales smaller than 10⁶, memory allocation dominates the time consumption, exceeding 90 % of the total execution time and significantly surpassing the durations of the other three processes. Notably, the time required for memory allocation remains approximately 0.4 s regardless of increases in data scale. The dominance of memory allocation at small data scales highlights its fixed overhead nature. This suggests that memory allocation is largely independent of data volume. While negligible in large-scale computations, this fixed cost becomes a critical bottleneck for small-scale tasks. When the data scale surpasses 10⁵, the overhead of memory copy rises rapidly, with a growth rate higher than those of kernel execution and synchronization processes. At a data scale of 10⁸, memory copy accounts for approximately 50 % of the total time and exhibits a tendency for further increase. Under this condition, the time contributions of memory allocation, kernel execution, and synchronization are approximately 14 %, 18 %, and 18 %, respectively. The rapid escalation of memory copy overhead underscores the limitations of host-device data transfer bandwidth. The growth rate of memory copy time implies that the data transfer between the CPU and GPU becomes one of the primary factors influencing the computational performance of programs in heterogeneous cluster environments, and the I/O-bound workloads often underutilize GPU compute capabilities. Future efforts could focus on reducing communication overhead through strategies such as unified memory architecture and asynchronous communication. Additionally, integrating mixed-precision methods (Váňa et al., 2017), converting variables with minimal impact on simulation results from double-precision to single-precision, could further enhance data transfer efficiency between CPUs and GPUs.

4.3.2 Coupling performance comparison

Following the integration of GPU-HADVPPM4HIP V1.0 into the EPICC-Model V1.6.0, computational efficiency on the EarthLab was improved through communication optimization described in Sect. 3.2 and enhanced thread and block collaborative indexing detailed in Sect. 3.3. Furthermore, the hybrid parallelization scheme outlined in Sect. 3.4 was employed to extend the parallel computing scalability of the EPICC-Model V1.6.0 on the EarthLab. As introduced in Sect. 3.1, three model versions were established: the baseline unoptimized version HIP-Ori, the communication-optimized version HIP-Opt1, and the collaboratively indexed version HIP-Opt2. Figure 12a displays the average elapsed time required for a 28 h simulation across these three versions. To ensure comparability, all tests adopted identical hardware configurations, including the MPI+HIP hybrid parallelization scheme with 10 CPU processes and 10 GPU-like accelerators, and were compiled using the −O3 optimization flag.

Figure 12 (a) The average hourly elapsed time required for a 28 h simulation across the HIP-Ori, HIP-Opt1, and HIP-Opt2 versions; (b) the average hourly elapsed time required for a 28 h simulation between the Fortran and HIP-Opt2 versions; (c) the hourly elapsed time required for a 28 h simulation across the original Fortran-based HADVPPM program on CPUs, the GPU-HADVPPM4HIP V1.0 with CPU-GPU data transfer, and the GPU-HADVPPM4HIP V1.0 without data transfer. The black dots represent the average values, and the unit is seconds (s).

The GPU-HADVPPM4HIP V1.0 module features triple nested loops, and each invocation by the EPICC-Model V1.6.0 necessitates data transfer between the CPU and GPU-like. Frequent CPU-GPU data transfer severely compromised the computational efficiency of the EPICC-Model V1.6.0 on the EarthLab. In the HIP-Ori unoptimized version, completing a 28 h simulation required 149.7 h, with an average hourly elapsed time of 19 241.0 s, equivalent to approximately 5.3 h, reflecting critically low efficiency. After implementing communication optimizations in the HIP-Opt1 version, the communication frequency between the CPU and GPU-like was reduced from roughly 9.8 million to 10. This optimization drastically decreased the total elapsed time from 149.7 h for HIP-Ori to 8.8 h for HIP-Opt1, while the average hourly elapsed time dropped from 19 241.0 to 1133.7 s, achieving a 17.0× improvement in computational efficiency. Subsequently, collaborative indexing of threads and blocks was applied to parallelize the two-dimensional grid computations in the EPICC-Model V1.6.0. Compared to HIP-Opt1, the HIP-Opt2 version further reduced the total elapsed time from 8.8 to 6.0 h, with the average hourly elapsed time decreasing from 1133.7 to 768.6 s. This enhancement delivered an additional 1.5× efficiency gain. Cumulatively, these optimizations elevated the computational efficiency of the EPICC-Model V1.6.0 on the EarthLab by approximately 25.0×.

By significantly improving the computational efficiency of the model on heterogeneous clusters through reducing the communication frequency between the CPU and GPU, we aim to emphasize the fundamental architectural differences between CPUs and GPUs. In recent years, although the substantial increase in GPU computational performance has opened up new directions for enhancing the efficiency of numerical models, the historical development of modeling frameworks means that early numerical models were predominantly designed for CPU architectures. Achieving parallel computation of numerical models on GPUs cannot be accomplished merely through straightforward code translation or by relying solely on heterogeneous programming interfaces. Instead, beyond adapting the model to GPUs, it is necessary to redesign the computational workflow, loop structures, and logical organization in accordance with GPU architectural characteristics, so as to fully leverage the powerful parallel computing capabilities of GPUs.

Figure 12b further compares the computational performance between the original Fortran-based version and the HIP-Opt2 version. In the figure, the label Fortran refers to the legacy CPU cluster implementation of the EPICC-Model V1.6.0, while HIP represents the HIP-Opt2 version. Following the parallelization of the advection module on GPU-like accelerators, the EPICC-Model V1.6.0 demonstrated superior computational efficiency on heterogeneous clusters compared to its Fortran counterpart on conventional CPU clusters. The total elapsed time for the 28 h simulation decreased from 22 874.0 s for the Fortran version to 21 520.0 s for the HIP-Opt2 version, while the average hourly elapsed time decreased from 817.0 to 768.6 s.

As demonstrated by the timing analysis of key processes in Sect. 4.3.1 for the offline heterogeneous computation of GPU-HADVPPM4HIP V1.0 on GPU-like accelerators, memory copying, specifically data transfer between the CPU and GPU, emerges as the dominant factor influencing parallel computational efficiency on heterogeneous clusters when handling large-scale datasets. This overhead surpasses the time spent on kernel function parallelization. To quantify this effect, the computational time of GPU-HADVPPM4HIP V1.0 was separately evaluated on GPU-like accelerators under two scenarios: (1) including CPU-GPU data transfer and (2) excluding data transfer (kernel-only computation). In Fig. 12c, the y-axis labels Fortran, HIP-IO, and HIP-without-IO correspond to the computational time of the original Fortran-based HADVPPM module on CPUs, the GPU-HADVPPM4HIP V1.0 with CPU-GPU data transfer on GPU-like accelerators, and the GPU-HADVPPM4HIP V1.0 without data transfer on GPU-like accelerators, respectively. All tests were compiled with the −O3 optimization flag, and timing metrics were averaged over a 28 h simulation, focusing on the hourly computational cost of the advection module.

The original Fortran-based advection module required an average of 110.6 s per simulated hour on CPUs. After heterogeneous porting and parallel optimization to GPU-like accelerators, the GPU-HADVPPM4HIP V1.0 with data transfer achieved an average time of 67.1 s per simulation hour, representing a 39.3 % improvement in computational efficiency. When excluding data transfer, the kernel-only GPU-HADVPPM4HIP V1.0 reduced the average time to 5.4 s per simulation hour, achieving a 20.5× acceleration compared to the Fortran version. These results underscore that while GPU-like accelerators deliver substantial computational power, the efficiency of CPU-GPU data transfer critically constrains overall performance on heterogeneous clusters. To mitigate this bottleneck, future work will focus on retaining I/O operations on CPUs while porting the entire physicochemical integration module (excluding I/O) to GPUs for parallel computation. This strategy is expected to reduce the impact of inter-device data transfer and further enhance scalability.

5 Conclusions and discussion

In recent years, the rapid advancement in the computational performance of GPUs has provided novel approaches and hardware foundations for improving the computational efficiency of air quality models. Building upon the heterogeneous porting and parallel optimization technology system for air quality model, this study further implements parallel computing of the advection module in the EPICC-Model air quality modeling system on GPU-like accelerators, validating the feasibility of the heterogeneous porting framework. The study involves three key technical improvements. The first is restructuring the original Fortran code of the advection module using standard C language programming. The second is porting the advection module to GPU-like accelerators through HIP heterogeneous programming technology, in addition, computational efficiency was enhanced through optimizing CPU-GPU data transfer frequency reduction, coordinated indexing of threads and blocks, and hybrid parallelization strategies. These optimizations collectively improved both the computational performance of the advection module and the parallel computing scalability of the EPICC-Model V1.6.0 on the EarthLab.

This study systematically evaluated the numerical consistency and computational efficiency of GPU-HADVPPM4HIP V1.0 through offline tests and further assessed the pollutant simulation performance of the EPICC-Model on EarthLab. Initial benchmarking compared the offline computational efficiency between the original Fortran-based HADVPPM program on domestic CPUs and the GPU-HADVPPM4HIP V1.0 implementation. The results demonstrated that the −O3 compiler optimization flag significantly enhanced GPU-HADVPPM4HIP's computational efficiency, with acceleration effects becoming more pronounced at larger data scales. Specifically, at 10⁸ data size configuration, GPU-HADVPPM4HIP V1.0 achieved a maximum 556.5× speedup over the Fortran baseline using default −O0 compilation, while maintaining a 128.0× speedup advantage even when both implementations employed −O3 optimization. Further profiling of GPU-HADVPPM4HIP's heterogeneous computation on GPU-like accelerators revealed critical characteristics: The elapsed time of memory-copy operations, namely CPU-GPU data transfers, increased rapidly with input-data size and accounted for approximately 50 % of the total computational time for an input array containing 10⁸ values. This observation underscores data transfer efficiency as a critical bottleneck for high-resolution air quality simulations in heterogeneous computing environments.

Coupling GPU-HADVPPM4HIP V1.0 into EPICC-Model V1.6.0 with data transfer optimizations and thread-block coordinated indexing strategies yielded system-level performance improvements of 17.0× and 1.5× respectively on the EarthLab. The detailed module-level analysis demonstrated that GPU-HADVPPM4HIP V1.0 achieved 39.3 % efficiency enhancement over the original Fortran advection module when accounting for CPU-GPU data transfer overheads, escalating to a 20.5× acceleration when excluding data transfer costs. These findings quantitatively validate the substantial impact of CPU-GPU data transfer efficiency on the operational performance of air quality models in heterogeneous computing architectures.

There remains considerable room for improving the computational performance of GPU-HADVPPM4HIP V1.0 within the EPICC-Model in this study. First, constrained by the thermal dissipation limits of transistors, the growth in computational performance of CPU processors has slowed and is gradually approaching its physical limits. In recent years, however, GPU processors have continued to achieve substantial performance gains due to their architectural advantages, and heterogeneous supercomputing architectures centered on “CPU+GPU” now dominate the landscape of advanced high-performance computing systems worldwide. Taking the 66th TOP500 list released in November 2025 as an example, 9 out of the top 10 supercomputers adopt heterogeneous architectures. Hence, heterogeneous computing represents the primary direction for the future development of supercomputing. That said, the advancement of numerical models relies fundamentally on the support of supercomputing capabilities. The team awarded the Gordon Bell Prize for Climate Modeling in November 2025 leveraged the powerful computational capacity of NVIDIA GH200 GPUs to accomplish the world's first 1.25 km ultra-high-resolution Earth system simulation (Klocke et al., 2025). To achieve kilometer-scale, ultra-high-resolution simulations of atmospheric pollution, it is both essential and imperative to adapt air quality models to heterogeneous supercomputing environments and enable efficient parallel computing.

Moreover, results from offline computational performance tests in this study and coupled performance tests in Cao et al. (2024) consistently show that the acceleration effect of GPUs becomes more pronounced as the computational scale increases, which sufficiently demonstrates the potential of GPUs in large-scale, high-resolution application scenarios. The simulation case configured in this study, with a resolution of only 15 km over the whole of China, may not yet fully reflect the parallel computing advantages of GPUs. In the future, higher-resolution simulation cases will be designed to better highlight the performance advantages of GPU acceleration.

Finally, in future work, we will employ further optimization techniques to enhance the computational performance of the advection module and other computationally intensive modules in the EPICC-Model on domestic heterogeneous clusters, including but not limited to:

1. First, future work will further improve CPU-GPU communication efficiency through chunked transfers using pinned memory, shared-memory and common-field caching, asynchronous communication, and unified memory. These strategies will reduce repeated data transfers and increase the overlap between CPU-GPU communication and GPU computation. In addition to these software-level optimizations, advances in hardware interconnects provide a complementary means of improving CPU-GPU communication efficiency. PCIe is a general-purpose, high-speed interconnect widely used for data exchange between host CPUs and accelerator devices. Its nominal per-lane raw data rate has increased from 16 GT s⁻¹ with PCIe 4.0 to 32 GT s⁻¹ with PCIe 5.0 and 64 GT s⁻¹ with PCIe 6.0, while the PCIe 7.0 specification further increases it to 128 GT s⁻¹ (Lender, 2019; Yanes, 2022; PCI-SIG, 2025). The EarthLab heterogeneous cluster used in this study is currently equipped with PCIe 4.0 buses; therefore, the performance improvement offered by PCIe 5.0 cannot be directly evaluated in the present study. If access to a heterogeneous computing platform equipped with PCIe 5.0 or a later-generation interconnect becomes available, we will evaluate GPU-HADVPPM4HIP V1.0 on that platform and further improve CPU-GPU communication efficiency through hardware-software co-optimization.

2. Second, while GPU-accelerated modules including the gas-phase chemistry module (Cao et al., 2025) and advection module have been individually developed, their systematic integration into EPICC-Model requires architectural refinement to increase GPU code coverage. We will analyze the computational characteristics of the code in other modules of the model. For code segments involving iterative computations, we will first decouple the iterative computations by creating intermediate variables, thereby eliminating dependencies between successive calculation steps. Subsequently, the computational code of other modules after iterative decoupling will be rewritten in the form of Kernel functions to increase the proportion of code executed on the GPU.

3. Finally, in current heterogeneous architecture supercomputing systems, the number of CPU processes within a computing node typically exceeds the number of GPUs. Employing the current matching scheme of one CPU process to one GPU accelerator results in the waste of remaining CPU computing resources. In the future, on one hand, while avoiding data transmission competition between the CPU and GPU, we will consider designing a more sophisticated mechanism for matching multiple CPU processes with a single GPU. On the other hand, drawing on Cao et al. (2024), we plan to introduce an OpenMP shared-memory parallel scheme into the EPICC-Model. Through multi-level hybrid parallelism, while porting computationally intensive modules to the GPU for parallel computing, other modules running on the CPU will utilize OpenMP multithreading parallelism, thereby fully leveraging CPU computing resources.

Consequently, achieving efficient parallel computation of air quality models on GPU accelerators necessitates collaboration between researchers well-versed in the intrinsic physicochemical mechanisms of the models and engineers with expertise in GPU hardware and software. Herein, we appeal for the participation of engineers specializing in GPU software and hardware development to jointly advance the progress of air quality modeling.