引言

数值预报作为气象预报的重要技术基础，它的发展与完善是提高气象预报整体水平的最有效途径。作为整个数值预报业务系统核心的中期业务数值预报，在中国的发展还远落后于欧洲中心、美国和加拿大等发达国家，主要表现在同化技术落后、模式分辨率低、物理过程的描述比较简单。因此，发展较先进的变分同化技术、提高模式的分辨率、细化和完善物理过程的描述，并采用先进高效的计算方法成为一项紧迫的任务。新一代中期数值预报模式以ECMWF的全球谱模式为基础，结合国家气象中心并行计算机的性能特点，以分布与共享相结合的方式实现了全球谱模式的高效运行。采用调整向量长度、优化程序设计、完善消息传递机制和实现MPI与OpenMP的混合并行编程等方法，减少了模式的通信量、计算量和内存的使用量，提高了计算效率，达到了业务对时限的要求。同时，设计开发了与模式相匹配的采用并行处理技术的全球资料分析同化系统和分析预报作业监控系统及资料预处理、模式后处理、产品生成、检验和绘图等外围系统。

本文主要介绍T213L31数值预报系统发展的过程中，所解决的主要技术问题、模式的特点和产品应用情况。

1 T213L31模式开发研制的技术难点 1.1 模式的移植和改造

引进的ECMWF业务版本的全球谱模式采用的是单任务多数据流的并行方式，并采用Fort ran 90和MPI消息传递库进行编程。由于此模式的编程是建立在富士通公司的VPP700向量计算机的基础上，而中国气象局1995年购置的IBM-SP2并行计算机是标量机，模式的移植主要要解决两个计算机之间的不兼容之处。另外，VPP700计算机采用向量处理机为计算节点，具有庞大的内存。为了使模式在VPP700计算机上高速运行，ECMWF采用了很长的向量长度。而我们的计算机是以基于超标量技术的工作站为基础建立的，长向量不仅降低了节点的计算速度，而且内存的占用量非常庞大，不作修改在SP2计算机上仅能计算T42L31以下的分辨率。因此，我们将长向量改造成适合于SP2的短向量，这不仅能够更好的发挥SP2高速缓存器的作用，使计算速度加快，同时还大大减少了模式内存的占用量，使在SP2上可以运行T106L31分辨率的模式。

1.2 模式的优化

1999年底中国气象局又购买了IBM SP机，其峰值计算性能达到800亿次/s，是32个处理器的IBM-SP2机 (85亿次/s) 峰值计算性能的近10倍，使运行T213L31分辨率的模式成为可能。但是，如果不做优化直接把模式从IBM-SP2机移植到SP计算机上，在T213L31分辨率下，虽然可以运行，运算速度比较慢，采用4个处理机进行10天预报需要42 h，无法满足业务对时限的要求。因此针对SP计算机的特点，对模式进行了以下几个方面的优化和改造:

(1) 计算编程的优化

SP计算机的处理机是基于高速缓存技术的超标量处理机，单处理机运算速度是系统整体性能的基础。每个处理机有5个浮点处理部件和两级高速缓存。提高高速缓存数据的利用率和充分发挥处理机内部5个浮点处理部件之间的并行性是提高运算速度的关键。我们采用的主要是CPU优化技术和内存优化技术。CPU优化技术包括回绕技术 (unrolling) 和多重加载技术 (multiple loads)，将循环内的计算操作分解成多个互不相干的计算操作，提高CPU内多个处理部件的并行处理效率；内存优化技术主要是减少对数据的存储操作，提高处理机内部的四个数据流水线的并行性，将数据分割成与高速缓存大小一致的数据块，使CPU反复使用数据块内的数据，从而增加对cache的命中率等。我们对模式中Legendre变换、短波辐射等计算量较大的12个子程序进行优化，使用了厂商提供的快速傅氏变换库函数，使模式的运算时间减少了一半，速度加快一倍。在4个处理机上进行10天预报时间缩短到20 h。原来在64个CPU上积分一个时间步长约需要18 s，经过优化后只需要9 s，可见优化效果是相当明显的。

(2) 消息传递机制的完善

消息传递是广泛应用于分布式可扩展并行计算机和工作站网络的编程模式，而MPI是实现这一编程模式的标准化的和可移植的消息传递库。原方案采用的是一种基于缓冲器的同步、阻塞式的消息传递机制，需要将数据在用户空间和系统缓冲区之间拷贝，增加了额外的内存之间的拷贝。因此，需要占用额外的系统内存空间，且增加了通信时间。同时，采用阻塞式的通信方式，数据拷贝到系统缓冲区 (对消息发送者) 或从系统缓冲区拷贝到用户缓冲区 (对消息接收者) 以后，处理机才能进行下一步的计算，这种方式虽然相对更安全，却增加了通信的时间。为了提高T213L31模式的加速比，减少通信所占的时间，我们将其通信方式改造成非同步、非阻塞的消息传递机制。采用这种方式，处理机只是启动一个发送或接收信息的操作，不用等通信完成，处理机即可返回继续下面的计算工作。因此使通信和计算时间重叠，从而隐藏通信延迟。

(3) 分布内存与共享内存并行的结合

SP计算机是由对称多处理节点组成的并行计算机系统，共有10个处理节点，每个节点由8个CPU组成。按照原方案，每一个CPU是一个单独的节点。但是在SP机上，如果使用快速通信协议 (US)，每个节点最多允许使用4个CPU，否则只能采用通信速度较慢的IP方式，因此必须使用共享内存并行与分布内存并行相结合的方式才能使用全部计算机资源。为此我们对模式程序进行改造，在节点内部使用OpenMP进行共享内存方式的多线程并行处理 (OpenMP是共享内存并行编程的工业标准)，在节点间仍然采用MPI消息传递的并行处理方式。实现了MPI与OpenMP的混合并行模式运算，大大减少了模式的通信量，减少了内存的使用量，提高了模式的运算速度。在64个CPU上T213L31模式积分一个时间步长由9 s减少到6 s左右。

经过上述开发，模式运算速度大大加快，在有限的计算节点下 (6个节点)，10天预报可以在2.5~3 h之内完成。同时，通过试验，我们发现采用OpenMP用2个线程，每个节点有4个MPI任务时，T213L31模式运算更快。即数据按尽可能相等的原则划分成为24个分区，分配给24个MPI任务 (每个节点4个)，这24个任务通过消息传递的方式交换信息；每个任务又在同一节点的2个CPU间采用共享内存的并行方式。图 1为在一个计算节点 (8个CPU) 上的任务并行方式。

1.3 并行化界面

原T213L31模式采用MPI作为消息传递的编程界面，但是模式并不直接调用MPI，而是通过一个中间的界面MPE调用MPI，如图 2所示。它的好处是很容易移植到使用其它消息传递界面的机器上，但缺点是增加了一层界面就增加了开销。在SP计算机上，为了采用非阻塞的消息通信机制，我们更改了部分程序，用模式直接调用MPI。

1.4 并行计算模式正确性检验

并行模式正确性的检验主要是通过比较用不同处理机数量时，模式计算结果有无变化，模式的各种统计量与在其它机器上得到的结果差异有多少，模式预报结果是否合理等方面来进行的。

(1) 不同处理机数量运行结果的比较在SP计算机上，分别用4，8，16，32，64个处理机进行12个时间步的预报，并对输出文件进行比较。我们发现二进制文件没有区别，说明预报结果完全一致。

(2) 统计量的比较将模式的1996年12月1日6 h预报的各层涡度、散度、动能、温度谱模的平均与ECMWF在VPP700上得到的相应结果进行对比，最大相对误差小于千分之一，可以认为预报结果是正确的。

(3) 预报结果合理性验证进行了1997年1月1日，4月1日，6月1日，7月1日，8月1日，10月1日的6天预报，并将预报结果与欧洲中心业务预报结果、实况和国家气象中心的业务预报进行对比。从预报结果可以看出，模式的预报结果与欧洲中心当时的T213L31业务预报比较相似。

经过以上检验，我们认为并行计算结果是正确的，T213L31模式已经成功地在IBMSP机上建立。

2 模式的基本特点 2.1 动力特点

(1) T213L31模式为三角截断的全球谱模式，截断波数为213个波，在格点空间其水平分辨率达到60 km，垂直方向为31个η面，模式顶为10 hPa，水平和对流层的垂直分辨率都比原来的业务模式T106L19提高了一倍。

(2) T213L31在格点空间采用规约化高斯格点，使计算量大大减少，内存需求和结果存储量也有比较大的减少。所谓规约化高斯格点是指减少赤道附近地区以外纬圈的高斯格点数，使这些纬圈的格点长度不超过在赤道的格点长度，同时剩下的格点数还能被用来进行快速傅氏变换。采用规约化高斯格点以后，覆盖全球的格点数能节省三分之一，模式计算精度并没有显著降低。

(3) T213L31采用半拉格朗日方法处理平流。欧拉方法是在空间中固定位置看流体的运动变化，拉格朗日平流的观点是跟着流体质点看流体运动。T106采用的是传统的欧拉平流方案，在此方法中，模式的积分时间步长受到平流风速的限制，在相同的计算稳定条件下，模式需要更多的积分步数。采用半拉格朗日方案后，积分时间步长不再受到平流风速的限制，积分时间步长可以增加5倍，T213L31的积分时间步长由欧拉方案的3 min增加到15 min。试验表明预报质量并没有显著的降低，而CPU时间减少了4倍。考虑到半拉格朗日方案的额外开销大约为20 %，采用半拉格朗日方案产生了相对于欧拉方案的4倍的有效率。

2.2 物理过程特点

与原业务模式T106L19相比，T213L31采用了一些更新、具有更真实物理概念的物理过程参数化方案，对于辐射方案、次网格尺度地形参数化、积云对流方案、云方案、陆面过程方案都作了较大改进。

T213模式的长波辐射是Morcret te^[1]的方案，短波辐射是Fouquart and Bonnel^[2]的方案；晴天长波通量的计算用比辐射率方法，同时用一个比T106L19更好的参数化方案来描述长波吸收对温度和压力的依赖关系。既考虑了水汽的p型连续吸收，又考虑了T106模式中没有的e型连续吸收。短波通量用光子路径分布方法，分开辐射传输中散射和吸收过程的贡献。散射的处理用Delta-Edding ton近似，透射函数用Pade近似。

T106采用了包络地形和Baines和Palmer^[3]的次网格尺度地形拖曳参数化方案 (重力波拖曳方案)；T213采用的是平均地形和一个新的次网格地形参数化方案 (Lott and Miller)^[4]。这个新的参数化方案既考虑了重力波拖曳，又考虑了阻塞流拖曳，与平均地形结合，克服了包络地形使同化过程中更多的低层数据被拒绝和人为提高地形所造成的过多对流降水的缺点。

在T106的积云对流参数化方案中，深对流采用Kuo^[5]的方案，浅对流采用垂直扩散方案 (Tiedtke)^[6]。T213采用的是Tiedtke^[7]的质量通量方案，引入了T106方案中没有的积云下沉支、积云动量传输和中层对流参数化，因此具有更真实的物理概念。

T106的云方案是Slingo^[8]的诊断云方案。T213采用的是Tiedtke^[9]的预报云方案。这个方案适当描写了次网格尺度凝结的动力影响，且云与辐射、动力和水文过程有更直接的联系。

T106的陆面过程方案把土壤分为三层，是Blondin^[10]的方案；这个方案基于两个活跃的土壤层的热量和水的收支方程，再加上一个土壤温度和湿度每月有个规定值的“气候层”，作为下边界条件。T213的陆表面参数化方案是Viterbo and Beljaar^[11]的方案，土壤被分为4层，4层的土壤温度和湿度都是预报量，热量和水分收支的下边界条件分别是零热通量和自由渗漏。它与T106的主要差别在于，考虑了周—季时间尺度的土壤水文过程，并且对土壤水文过程进行了更物理化的描述。

3 客观分析方案的开发 3.1 T213观测资料的生成

建立了独立于Cray-C92计算机的两套观测资料预处理系统，它们分别与Digital Alpha机的要素库系统和9210商用数据库相连。对两套系统中观测资料的数量、类型及对分析的影响等做了大量的对比试验工作，两者差别不大。目前业务采用要素库系统为T213分析程序提供观测资料。

3.2 全球客观分析方案的开发

全球客观分析方案是全球中期数值预报业务系统的一个重要组成部分，它将全球可以使用的观测资料进行分析同化处理，为全球业务预报系统提供数值形式的初值。T213L31模式引进时没有引进配套的客观分析方案，因此需设计开发与T213L31模式匹配的分析方案。T213L31分析方案的开发建立在原T106L19最优插值 (OI) 分析方案的基础上，包括移植、升级和并行计算三个主要方面。

首先，将原T106L19最优插值 (OI) 串行分析程序从CRAY-C92机上移植到SP机上，去除原有的CRAY机器上使用的POIN TER语句等非标准语句，修改了原有的内存管理系统，补充所缺子程序 (原调用CRAY子程序库)，去掉程序中的工作文件，以加快运算速度。然后，对移植后的T106L19 OI串行分析程序进行升级，使其与模式分辨率相匹配，即对水平方向213波三角截断的谱系数、垂直方向的31个模式层进行分析。最后，对升级的分析程序进行并行计算开发，主要是针对耗时大的子程序 (占总计算量的90 %) 进行并行化处理。采用了MPI消息传递库，保证了程序的可移植性。

开发客观分析方案的并行程序的工作大致分为三步:确定客观分析方案的并行计算方案；开发并调试并行程序；进行试算和正确性验证，分析其运行性能并加以优化。

客观分析程序将全球分成许多小的“方盒子”，每一个“方盒子”之间的计算是相互独立的，每个“盒子”的计算量取决于盒子内观测资料的多少。但是由于观测资料在空间上分布极不均匀，而且不同时间的资料情况也不相同 (图略)。在无法事先知道计算量的情况下，采用固定分配“盒子”数量的方法会造成节点间计算负载的不平衡。我们采用了主—从方式的并行处理方法，主节点只负责将工作分配给从节点，从节点接受到任务后开始计算，计算完成后将结果传回主节点，主节点接到结果后，再传给它下一个任务直到完成。这样，尽管每个从节点的速度有差异，工作量可不相同，但都处在忙的状态，实现了工作负载的最佳平衡。图 3是具体的流程图，其中箭头表示两者间存在信息通讯。MAS TER按原串行程序流程执行，但在运行到资料检查、分析系数估算、格点估算三个子程序时，只负责控制任务的分配。S LAVE具体执行这三个子程，但其数据来源于MASTER。

由于并行的三段子程序在串行运算过程中所占时间为90 %，而且三段子程序的SLAVE并行全部具有自动负载平衡功能，因此并行程序的效率是非常高的。用若干时次资料进行试算，其并行程序与串行程序的结果差异很小。因此，我们认为分析并行程序开发是成功的。

另外，经过对2个CPU、4个CPU和8个CPU的并行运行对比，综合考虑程序的并行效率、机时条件和任务要求，我们认为对于T213的分辨率，使用8个CPU并行计算比较合适。图 4给出的是全球分析方案在SP2机上的运行情况。图中所标的并行和串行的数字，分别是模式的并行运算和串行运算所用的时间。

3.3 矩阵病态现象的解决

在T213L31分析方案试运行中，我们发现有时会出现矩阵病态现象，使分析无法进行下去，因此不能为模式提供初值。对观测资料分析后发现，产生问题的原因是个别观测资料之间的高相关性。因此从数学角度采取了调整线性方程组重新求解，即矩阵降阶技术处理 (在物理上即是去除个别资料，对实际应用效果影响不大，可认为是合理的)，以保证方案运行的可靠性和稳定性。

4 预报同化系统的建立

预报同化系统是指每天4次的间歇性6 h同化过程和模式在00 :00 UTC的3天和12 :00 UTC的10天预报。这个系统流程的开发是由我们自己独立完成的，实现了最优插值方案与T213L31模式的集成。由于模式的全部输入输出都采用GRIB码文件，好处是节省空间而且其内部含有对数据的描述，但是分析程序需要IEEE文件。因此我们开发了OI分析后高度场、湿度场、涡度、散度和地面气压场压缩成GRIB码程序，模式气候资料场生成程序，和由模式输出的GRIB码解压缩到OI分析程序需要的高度、涡度、散度和地面气压场的数据接口处理程序。另外，原业务计算机平台CRAY-C92是只有两个CPU的共享并行矢量机，IBM SP是大规模并行的分布与共享标量机，在文件系统、作业控制和管理方面与前者都有很大的不同。因此，我们研制开发了适合大规模并行计算机的，由分析到预报再到分析的自动化作业流程，对机器资源的申请、临时文件的管理、运行监控、异常情况的处理等都进行了考虑，解决了4次连续分析同化预报的自动化流程控制和运行监控的技术问题。图 5是T213L31模式的同化预报流程图。

5 后处理系统的自行开发

T213L31模式引进时没有后处理系统，因此必须自行研制开发能满足业务需求的后处理系统。图 6是我们开发的后处理系统的流程图。

在试验的过程中我们还发现，由于T213L31的资料量远比T106L19多，串行的后处理需要6、7 h才能完成，无法满足业务对时间的限制。针对这个问题，对后处理部分进行了多任务并行处理，完成这部分工作的墙钟时间减少到不到15 min，这样不但能缩短整个系统的运行时间，尽早为预报员提供产品，还充分利用了计算机资源。

6 产品的应用情况

自从T213中期预报系统试运行以来，国家气象中心数值室一直对模式预报进行跟踪检验和效果评估。这期间的统计学的检验表明，无论是北半球还是东亚地区的形势场检验，T213各时效的500 hPa距平相关系数基本上都高于相同时效的T106预报；在冬、春季，T213的北半球预报有效时效比T106延长了1天多，即使是在夏季的6月份，预报有效时效也延长了近1天；对于东亚区来说，也有类似的特点。T213对降水落区 (≥0.1mm) 的预报比T106有了明显的改善，空报率比T106减少。对较大量级降水的预报，Ts评分也比T106高。天气学的检验表明，T213对西风指数的逐日演变、副高脊点和北界的进退变化和48 h之内台风的移动变化预报都较T106有了明显的改善；T213模式对降水落区或雨带位置的预报优于T106和HLAFS，3天以内的预报有较高的可信度。对于其产品进行试用的预报员和其他数值预报产品的使用者也都反映这个新业务系统的预报产品比原T106系统有了明显的改进。

由于篇幅的原因，这里只给出2001年中国北方最大的一次降雨过程，T213系统的预报情况。图 7给出的是T213对这次过程连续两天的24 h累积降水预报和相应时段的实况降水分布。可以看出无论是对降水落区还是降水中心的预报，T213的预报都与实况有较好的吻合。

7 结束语

T213L31业务中期数值预报系统的建立缩小了我国的中期数值预报与发达国家的差距，为预报员的预报提供了更好的参考。但同时，我们也认识到进一步发展我国的中期数值预报业务的迫切需求。目前，欧洲中心正在开发T799L90模式，很快将会实现业务化。日本气象厅的全球模式2006年将由目前的业务模式T213L40，升级到T959L60。加拿大气象局的全球模式将由100 km ×100 km-28层升级到35 km ×35 km-60层。因此，进一步提高我国中期数值预报模式的分辨率和进行相应的物理过程的改进工作，是必然的也是很紧迫的一项任务。国家气象中心2003~2005年事业发展计划也对此提出了要求，目标是2005年建立Tl319L31全球中期数值天气预报系统。另外，这个发展计划指出建立水平分辨率为几百米的奥运精细数值预报模式，该模式的建立也需要一个比现有业务系统更高分辨的全球模式提供边值条件。

致谢 本文是对许多同志共同开发结果的总结和介绍。T213L31中期业务预报系统在IBM SP机上的建立工作在陈德辉、王建捷、石曙卫、张跃堂的先后组织领导下，国家气象中心数值室其它同志如，张德新、万丰和赵刚对客观分析系统的开发，胡江凯对作业流程和运行监控的建立，赵俊英、王雨和黄卓对检验系统的建立，庄建敏和苏颖对场库的建立和完善，王毅涛对micaps绘图的研制，管成功、马秀君、陆志善、王克敏和李纪曼对与系统业务化相关的工作都做出了各自的贡献，也为本文的完成提供了一些资料。本文的完成还要感谢陈德辉博士，文中关于加拿大和日本气象厅的发展计划来自于其参加东京暴雨模式国际报告会后的介绍。另外，在此还要感谢IBM公司的Zaphiris Christidis先生对T213L31中期业务预报系统建立工作的帮助和支持。