# Using a GPU to Accelerate a Longwave Radiative Transfer Model with Efficient CUDA-Based Methods
## 4 GPU-Enabled Acceleration Algorithms
### Acceleration Algorithm with Two-Dimensional Domain Decomposition  
\(nlay\), the number of model layers in the vertical direction of the RRTMG_LW, is 51 in the CAS-ESM. Besides the parallelization in the horizontal direction, the one in the vertical direction for the RRTMG_LW is also considered to make full use of the performance of the GPU. Figure 4 illustrates the two-dimensional (2-D) domain decomposition for the RRTMG_LW accelerated on the GPU. The acceleration algorithm in the horizontal and vertical directions or with a 2-D domain decomposition for the RRTMG_LW is illustrated in Algorithm 4. Here, \(tBlock\) defines the number of threads used in each thread block of the \(x\), \(y\), and \(z\) dimensions by the derived type \(dim3\). \(grid\) defines the number of blocks in the \(x\), \(y\), and \(z\) dimensions by \(dim3\). Because of data dependency, \(cldpmc\) and \(rtrnmc\) are unsuitable for the parallelization in the vertical direction. These kernels with 2-D decomposition in Algorithm 4 are described as follows:  
(1) The implementation of CUDA-based \(inatm\_d\) with 2-D decomposition is illustrated in Algorithm 5. Here, \(threadIdx\%x\) identifies a unique thread inside a thread block in the \(x\) dimension, \(blockIdx\%x\) identifies a unique thread block inside a kernel grid in the \(x\) dimension, and \(blockDim\%x\) identifies the number of threads in a thread block in the \(x\) dimension. In the same way, \(threadIdx\%y\), \(blockIdx\%y\), and \(blockDim\%y\) identify the corresponding configuration in the \(y\) dimension. \(lay\), the coordinate of the vertical direction, also represents the ID of a global thread in the \(y\) dimension, which can be expressed as \(lay=(blockIdx\%y-1)\times blockDim\%y+threadIdx\%y\). Therefore, a grid point in the horizontal and vertical directions can be identified by the \(iplon\) and \(lay\) variables as two-dimensional linear data. Due to data dependency, a piece of code in \(inatm\_d\) can be parallel only in the horizontal direction, so the kernel \(inatm\_d3\) in Algorithm 4 uses 1-D decomposition. The other code in \(inatm\_d\) is able to use 2-D decomposition. Due to the requirement of data synchronization, \(inatm\_d\) with the 2-D decomposition is divided into four kernels (\(inatm\_d1\), \(inatm\_d2\), \(inatm\_d3\), and \(inatm\_d4\)). The 2-D parallelization of \(inatm\_d1\), \(inatm\_d2\), and \(inatm\_d4\) is similar to Algorithm 5. Their detailed implementations will not be described further here. According to the design shown in Algorithm 5, \(ncol\times nlayers\) threads can execute 2-D \(inatm\_d\) concurrently. (2) Similarly, due to data dependency, a piece of code in \(setcoef\_d\) can be parallel only in the horizontal direction, so the kernel \(setcoef\_d2\) in Algorithm 4 uses 1-D decomposition. The other code in \(setcoef\_d\) is able to use 2-D decomposition, as shown in the kernel \(setcoef\_d1\). The 2-D parallelization of \(setcoef\_d1\) is described in rough form in Algorithm A5 of Appendix A.2. (3) The 2-D acceleration algorithm of \(taumol\_d\) is illustrated in Algorithm A6 of Appendix A.2. Here, the 2-D parallelization of \(taugb1\_d\) with the device attribute is also described. The 2-D parallelization of the other 15 subroutines (\(taugb2\_d\small{\sim}taugb16\_d\)) is similar to \(taugb1\_d\). Figure 4: Schematic diagram of 2-D decomposition for the RRTMG_LW in the GPU acceleration. ```
subroutine_rrtmg_lw_d2(parameters)
1.