引言

随着计算机和计算方法的飞速发展，大部分学科通过采用定量化的方法变得更加精确，计算机成为科学研究等不可缺少的工具。与此同时人们对计算机的性能要求越来越高，而任何高性能计算和超级计算都离不开使用并行技术^[1]。

社会经济的发展对气象预报的要求越来越高，公众需要定点、定时、定量的气象服务，高分辨率数值预报模式是提供这种需求有效的方法之一。由于高分辨率数值预报模式计算量非常大，同时要求预报的时效性，对计算机性能要求非常高，在并行计算机上采用并行计算技术可以有效地解决数值预报模式计算量大、时效性强这个问题。

广州区域数值预报模式是广州热带海洋气象研究所自主开发的国内最早用于业务预报的有限区域原始方程模式^[2]，也是中国台风专家组准许发布台风预报的数值模式之一。早期的模式分辨率较低，物理过程也不完善。为了满足业务需要，模式在空间分辨率和物理过程等方面进行了升级，升级后模式垂直方向分为10层，水平分辨率为111.1 km (~1°)，时间积分采用省时显式积分方案 (EES)^[3]，对热带气旋中心位置的预报有了明显的改进^[4]。为进一步提高模式的业务预报能力，特别是提高气旋路径及强度的预报，最近又提高了模式的空间分辨率，使模式垂直方向达到22层，水平分辨率达到55.6 km (~0.5°)；同时模式预报区域也要扩大，以适应预报时间由原来的48 h延长到72 h，减少模式边界的影响。这使得模式的计算量迅速增加，串行模式已经不能满足业务的时效性要求，因此模式必须采取并行化计算。计算机条件的改善也为并行化的实现提供了基础。本文针对广州区域数值预报模式的并行化计算进行了一些研究和探讨。

1.   模式基本介绍

1.1   模式基本方程

模式采用σ坐标，水平方向采用Mercator投影，模式的水平动量方程为

模式热力学方程为

模式连续方程为

模式状态方程为

其中:

1.2   时间积分方案

模式起步采用Euler后差积分，即

这里f^{*n +1}是采用向前差分求得 (n +1) 时刻的第一近似值，f^{n +1}是采用向后差分求得 (n +1) 时刻的值，是指由第一近似值f^{*n +1}所计算的倾向值。

模式大部分时间积分采用省时显式积分方案 (EES)。设定两个时间步长Δt_a和Δt_b，Δt_a=m·Δt_b，m取4或6，时间积分按Δt_b进行，但所有的非线性项在2Δt_a时段内保持不变。以τ表示短步长的时间变量，向量x=(u，v，t，π，Λ) 表示预报量值，则有

这里N表示非线性项，而L表示线性项。长时步与短时步预报结果的联结关系是

这里xⁿ(m) 表示短时步第m个时间层上的值。

短时步积分过程中，速度和高度、温度定义在不同的时间层上，并采用“时间跳点”方案。这时有

其中

预报时和项在整个长时步2Δt_a内保持不变。

1.3   模式垂直网格和水平网格

模式采用σ坐标，如图 1所示，实线代表整层，虚线代表半层。模式水平方向采用半跳点的Arakawa B网格，如图 2所示，图中V表示矢量格点如u，v；F表示标量格点如，T′，π，。

图  1  模式的垂直分层示意图

下载: 全尺寸图片 幻灯片

图  2  模式变量水平网格示意图

下载: 全尺寸图片 幻灯片

1.4   模式侧边界

模式侧边界处理采用开放边界方法和边界松弛方法。

开放边界方法，即在边界上假设

其中C_x、C_y根据内点的变化倾向来估计。

边界松弛方法即在模式边界采用如下处理

这里是变量f的最终预报倾向值，表示更大尺度模式提供的预报倾向值，表示模式预报的倾向值，α是边界平滑权重系数。

2.   并行算法设计

模式的时间积分均为显式积分，即变量后一时步的求解只利用到前面时步的已知值，因此模式采用并行化计算从理论上可行^[5]。模式侧边界处理也采用的是显式方案，同样可以采用并行化计算，侧边界处理采用先进行数据传递后进行计算的方法。

模式在垂直各层上计算与水平方向上的时间积分不相关，因此模式可以在水平方向上采用分区并行计算，而在垂直方向上不适合进行分区并行。模式中最高偏导数为二阶偏导，采用格点差分求解时，某一内部格点上值的求解一般只用到周围8个格点上的值^[6]，在分区时区域内部的分区边界至少要外扩一圈，模式中的标量和矢量采用两套不同的网格。图 3(a)和3(b)分别为模式一维Y方向分区和二维X-Y方向分区并行方案。

图  3  模式一维Y方向分区并行方案 (a) 和二维X-Y方向分区并行方案 (b)

下载: 全尺寸图片 幻灯片

如图 3(a)和3(b)所示，双箭头表示每次积分时各分区边界要进行数据的相互交换情况。模式一维Y方向分区方案在并行方法上相对比较简单，数据交换时仅交换边界上的一行数据。这里一维并行设计采用自动定义各分区长度的方法，并行计算时可以任选一个到多个CPU运行，但要求每个分区在X和Y方向均不少于5个格点。一般情况下模式在水平方向上的格点数是固定的，采用一维分区并行计算时各分区在Y方向上的格点长度不一定相同，有的分区可能会多一行格点。

二维X-Y方向分区方案跟一维Y方向分区方案基本相同，只是它在X和Y方向同时进行分区。采用二维分区并行计算时各分区在X和Y方向上的格点长度不一定相同，有的分区可能多一行或多一列格点。采用二维分区方案在数据交换时比一维分区方案复杂很多，其中有的分区要进行四个面的边界数据交换，因此模式实现二维并行比实现一维并行要复杂。另外，本模式采用二维分区并行方案时，通讯时间占总运算时间的百分比可能会随着分区数增加而明显增加。

模式采用一维Y方向分区和二维X-Y方向分区并行计算时模式积分和计算过程中没有采用隐式和迭代计算，模式并行计算只进行了分区处理，因此模式并行计算和串行计算的结果应该完全相同。这里把模式并行计算和串行计算的所有输出预报结果采用二进制文件进行比较，结果显示模式并行和串行计算输出的二进制文件完全相同，这表明模式并行和串行计算结果完全相同。这里模式并行采用MPI (Message Passing Interface) 并行计算环境。

3.   试验结果及并行效率分析

这里测试的并行模式格点数为:181 ×109 ×22，水平范围为:75.0°~165.0°E和0.0°~47.4°N，水平分辨率为55.6 km (~0.5°)。模式积分步长为120 s，积分时间为72 h。这里所有测试采用同一个例。

试验采用广州曙光3000并行计算机，该并行计算机拥有4个节点，每个节点有4个CPU共享内存，每个CPU为1.5 Gflops、2.0 G内存，节点间通讯网为1000 M/s。模式一维Y方向分区和二维X-Y方向分区的并行测试结果如表 1所示。

表  1  广州曙光3000并行机上模式一维Y方向分区和二维X-Y方向分区的并行测试结果

下载: 导出CSV 

| 显示表格

这里，加速比s_p的定义为

其中，T_l是单机上求解一个问题所需要的时间；T_p是在具有p个 (与单机相同) 处理机的并行系统上，并行求解同一问题所需要的时间。并行效率E_p的定义为

其中，s_p为p个处理机的加速比；p为处理机的个数。

通讯时间百分比的定义为

这里，总运算时间由通讯时间和实际计算时间两部分组成，其中通讯时间通过并行计算中的MPI-WTIME () 计时函数累计模式中数据通讯所用的时间得到。

表 1给出了广州曙光3000并行计算机上广州区域数值预报模式两种并行化方案的测试结果。可以看出:模式并行化以后采用多个CPU并行计算比模式串行计算的计算时间明显缩短；模式在使用相同个数的CPU时，一维分区并行方案比二维分区并行方案的并行效率要高。另外，模式采用一维分区并行方案时使用2~8个CPU的并行效率比较高，而使用12~16个CPU时由于通讯时间百分比的增加，模式并行效率明显降低，其主要是由于使用12~16个CPU时通讯时间百分比超过了50 %，即通讯时间比实际计算时间要长，模式并行运算的时间主要花在了数据通讯上面。这里值得注意的是模式采用二维分区并行方案时，使用12~16个CPU时比使用8个CPU时的运算时间要长，而且使用CPU的个数越多运算时间越长，出现这种现象的直接原因是由于CPU个数的增加，通讯时间明显增加，导致了总运算时间的增加，这可能跟模式的二维分区并行技术方案等有关。从并行测试结果可以看到，模式并行时不一定是使用CPU个数越多越好，只有在通讯时间百分比比较小的情况下，模式的并行效果才比较好。

广州区域数值预报模式采用一维分区并行化方案相对比较简单，目前该并行方案已经在数值预报业务中得到应用。一般模式并行效率跟模式本身的特点和模式的并行技术方案等有关，因此如何进一步减少模式通讯时间百分比提高模式并行效率以及改进模式二维分区并行方案将有待于进一步探讨。

4.   结论

(1) 广州区域数值预报模式使用的是显式积分和计算，并且模式在水平方向上的计算和垂直方向上的计算相互独立，因此模式适合采用水平分区并行计算方法。模式并行化以后采用多个CPU并行计算比模式串行计算的计算时间明显缩短，模式采用水平方向一维分区并行方案相对比较简单，而采用二维分区并行方案比较复杂。

(2) 模式采用水平方向分区并行计算时在模式积分和计算过程中没有使用水平方向上的隐式或迭代求解，所以模式并行计算结果和串行计算结果完全相同。模式采用一维分区方案比二维分区方案并行效率要高，模式两种并行方案在使用8个以上CPU时由于通讯百分比明显增加，并行效率显著下降。一般模式并行效率跟模式本身的特点和模式的并行技术方案等有关，在模式并行时要提高模式并行效率，就要尽量减少通讯时间百分比。模式采用一维分区并行方案在使用8个CPU时并行效率保持在90 %左右，可以满足业务需要。