# Using a GPU to Accelerate a Longwave Radiative Transfer Model with Efficient CUDA-Based Methods
## 4 GPU-Enabled Acceleration Algorithms
### Acceleration Algorithm with One-Dimensional Domain Decomposition  
The computation of the RRTMG_LW in the horizontal direction is independent, so each CUDA thread is assigned the computation of one column. It means that the RRTMG_LW is able to be parallel in the horizontal direction. It is noteworthy that the first dimension of 3-D arrays in the CAS-ESM RRTMG_LW CUDA code represents the number of horizontal columns, the second dimension represents the number of model layers, and the third dimension represents the number of \(g\) points. Figure 3 illustrates the one-dimensional (1-D) domain decomposition for the RRTMG_LW accelerated on the GPU. The acceleration algorithm in the horizontal direction or with a 1-D decomposition is illustrated in Algorithm 2. Here, \(n\) is the number of threads in each thread block and \(m=\lceil(\textbf{real})ncol/n\rceil\) is the number of blocks used in each kernel grid. In Algorithm 1, \(inatm\), \(cldprmc\), \(setcoef\), \(taumol\), and \(rtrnmc\) are all called repeatedly for \(ncol\) times. However, due to using \(ncol\) CUDA threads to execute these subroutines, they are only called once in Algorithm 2. Theoretically, their computing time will be reduced by \(ncol\) times. These kernels in Algorithm 2 are described as follows:  
1. The implementation of the kernel \(inatm\_d\) based on CUDA Fortran is illustrated in Algorithm 3. Within, \(threadIdx\%x\) identifies a unique thread inside a thread block, \(blockIdx\%x\) identifies a unique thread block inside a kernel grid, and \(blockDim\%x\) identifies the number of threads in a thread block. \(iplon\), the coordinate of the horizontal grid points, also represents the ID of a global thread in CUDA that can be expressed as \(iplon=(blockIdx\%x-1)\times blockDim\%x+threadIdx\%x\). According to the design shown in Algorithm 3, \(ncol\) threads can execute \(inatm\_d\) concurrently. 2. The other kernels are designed in the same way, as shown in Algorithms A1-A4 of Appendix A.1. To avoid needless repetition, these algorithms are described only in rough form. 3. \(taumol\_d\) needs to call 16 subroutines (\(taugb1\_d\small{\sim}taugb16\_d\)). In Algorithm A3 of Appendix A.1, \(taugb1\_d\) with the device attribute is described. With respect to the algorithm, the other 15 subroutines (\(taugb2\_d\small{\sim}taugb16\_d\)) are similar to \(taugb1\_d\), so they are not described further here. ```
subroutine\(rtmg\_lw\_d1\)(parameters)
1. Copy input data to GPU device //Call the kernel \(inatm\_d\)
2. call\(inatm\_d\lll m,n\gg\)(parameters) //Call the kernel \(cldprmc\_d\)
3. call\(cldprmc\_d\lll m,n\gg\)(parameters) //Call the kernel \(setcoef\_d\)
4. call\(setcoef\_d\ll m,n\gg\)(parameters) //Call the kernel \(taumol\_d\)
5. call\(taumol\_d\lll m,n\gg\)(parameters) //Call the kernel \(rtrnmc\_d\)
6. call\(rtrnmc\_d\lll m,n\gg\)(parameters)
7. Copy result to host //Judge whether atmospheric horizontal profile data is completed
8. if it is not completed goto 1 end subroutine
```  
**Algorithm 2**CUDA-based acceleration algorithm of RRTMG_LW with 1-D domain decomposition and the implementation of 1-D GPU version of the RRTMG_LW (G-RRTMG_LW). ```
**Algorithm 3**Implementation of \(\mathit{inatm\_d}\) based on CUDA Fortran in 1-D decomposition. ```  
attributes(global) subroutine(\(\mathit{inatm\_d}\)(parameters)  
1. \(\mathrm{iplon}\)=(\(\mathrm{blockIdx\%x\_1}\)) \(\times\)\(\mathrm{blockDim\%x\_4}\)
2. \(\mathbf{if}\)(\(\mathrm{iplon}\)\(\geq\)\(1\).and.\(\mathrm{iplon}\)\(\leq\)\(\mathrm{ncol}\))
3. Initialize variable arrays //nlayers=nlay + 1, nlay is number of model layers
4. \(\mathbf{do}\) lay = 1, nlayers-1
5. \(\mathrm{paved}\)(\(\mathrm{iplon}\),lay) = \(\mathrm{play\_d}\)(\(\mathrm{iplon}\),nlayers-lay)
6. \(\mathrm{tavel}\)(\(\mathrm{iplon}\),lay) = \(\mathrm{play\_d}\)(\(\mathrm{iplon}\),nlayers-lay)
7. \(\mathrm{pz}\)(\(\mathrm{iplon}\),lay) = \(\mathrm{plev\_d}\)(\(\mathrm{iplon}\),nlayers-lay)
8. \(\mathrm{tz}\)(\(\mathrm{iplon}\),lay) = \(\mathrm{tlev\_d}\)(\(\mathrm{iplon}\),nlayers-lay)
9.