《1 前言》

1 前言

气象灾害（ 水灾、旱灾、风灾、雪灾等）与地表灾害（ 海啸、泥石流、潮灾等） 、地质构造灾害（地震、火山爆发等）、生物灾害（ 农作物和森林的病虫害、鼠害、传染病流行等） 等其他类型的自然灾害相比，无论在影响范围的广度与深度、影响社会的人数与活动等方面，还是在造成直接（ 间接）经济损失方面，都是最大的，而且一年四季都可能发生。 随着社会经济的发展， 气象灾 害造成 的直接 经济损 失也在加大。

准确的天气预报是积极防御气象灾害、减少人员和财产损失、间接提高生产力的重要手段，而现代天气预报则是以数值天气预报技术作为重要科技支撑。 所谓数值天气预报，通俗地说就是计算机天气预报。 也就是说，根据大气演变的物理规律（如人们熟知的牛顿运动定律、热力学定律、质量守恒定律等），应用数值计算的方法，建立数学计算模型，编写成计算机能自动运行计算的软件，安装在巨型计算机上，每天定时输入气象观测数据，计算机就可以在很短的时间内制作出未来几小时、几天、任意地点的天气预报。 因此，建立在现代探测技 术、通信技术、计算机技术以及大气科学理论基础上的数值天气预报是解决天气预报准确率持续提高，满足社会经济发展对防灾减灾气象保障服务需求的根本科学途径。 数值天气预报的思想 最早是由挪威科学家Bjecknes 于 1904 提出来的[1] ，但是直至第二次世界大战结束后，真正意义的“计算机预报天气” 才获得成功[2]  。 自 此以后，数值天气预报准确率不断提高，在现代天气预报业务中发挥着越来越重要的支撑作用。

以世界上最先进的欧洲中期预报中心（ECMWF）为例，来说明业务数值天气预报水平的变化（见表 1）。从表 1 可以看出，从 1980 年至 2010 年，500 hPa 高度场预报距平相关系数*倡逐年都在提高， 比如预报时效为 3 d 的距平相关系数 1980 年为 84 %左右（北半球），2010 年已高达 98 %，第 5 天、第 7天、第 10 天的距平相关系数也同样得到提高， 这说明数值天气预报技术的科学性及其支撑能力在不断提升。

《表1》

表1 ECMWF 全球模式的北半球 500 hPa 高度场预报距平相关系数

Table 1 Anomaly correlation coefficient at 500 hPa height field for northern hemisphere at ECMWF

注：数据根据 Ro bert o Buiza 博士 2011 年 11 月 21—25 日访问中国气象局数值预报中心作的学术报告整理而得

正是基于数值天气预报技术在现代天气预报业务中的重要科技支撑作用，世界各国都投入大量的人力物力，布设地面、雷达、高空、卫星等气象探测网，购买巨型高性能计算机，发展先进的数值天气预报系统技术，为气象灾害天气预报预警提供有力支撑。中国气象局在国家科技研究计划的支持下，自 2001 年起，自主研制发展了我国的新一代全球与区域一体化数值天气预报系统（GRAPES）[3,4] 。

《2 GRAPES 模式的科学设计策略》

2 GRAPES 模式的科学设计策略

地球表面的大气运动变化（ 也即天气变化） ，遵循基本的物理定律：牛顿运动定律、热力学定律、质量守恒定律和大气状态变化关系[3～5]。

大气运动方程（牛顿运动定律） ：

式(1) ~(3) 中， 分别为东西、南北、垂直 风速分量； r 为球半径（可近似取为地球半径） ； 和 为科氏力项； λ 为经度； φ为纬度； t 为时间；p 为气压； δ1 和 δ2 分别为“ 0 ”或“ 1” 参数； Fx 如 Fu 等为次网格物理过程参数化项。

连续方程（质量守恒定律）：

式（4）中，为净热源热汇项；= 为三维散度；（为定压比热， R 为干空气体常数）； θ 为位温； 为气压变量。

热力学方程（热力学定律）：

描述温度 T 、气压 p 和密度 ρ 之间变化关系的大气状态方程为：

式（1） ～（8）即为方程组。 自 20世纪 50 年代中期以来的 60 多年进程中[6]， 以方程组为基础，各气象业务中心根据不同的预报对象需要，研究发展了多套不同的数值天气预报模式，并同时运行中尺度模式、区域暴雨模式、台风模式、全球预报模式、短期气候预测模式、沙尘暴模式等。 这些模式的差异性很大，使得系统运行维护与改进开发的成本高，研究成果的业务转化周期长[7,8]。 近年来世界各国着手研究建立区域模式与全球模式统一、天气与气候（ 月、季、年） 一体化模式、研究模式与业务模式一体化的新一代多尺度统一数值预报模式，以取代现有的多套不同尺度、不同预报对象的模式共存的数值预报业务体系，大大降低运行成本，加速业务成果的应用转化。 目前，加拿大[7,8] 、英国[9] 、法国[10] 等国家的气象业务中心都建立了程度不同的多尺度统一模式系统。

自 2001 年起，中国气象局根据我国气象数值天气预报业务发展的实际需要，持续实施了研究发展我国新一代全球与区域一体化数值天气预报系统 GRAPES 的科技攻关／支撑计划。

《2.1 模式大气近似假设》

2.1 模式大气近似假设

式( 1 ) ~( 6 )方程组是一组不做任何简化的“完整”的 方程组，它几乎包含大气中各种尺度的运动，同时考虑了三维科氏力、球面半径与球面曲率等影响。 GRAPES 模式采用浅薄的（ r ≡ =6 371 km） 、“ 全可压的” （ δ1 =1 ）非静力平衡大 气的近似假设（见表 2 ） 。

《表2》

表2 大气模式方程组的不同简化假设方案

Table 2 The simple hypothesis methods for equations of the atmospheric model

注：将 r ≡ =6 371 km 时的模式大气定义为“浅薄大气” ，将 r > 时的模式大气定义为“ 深厚大气” ，并给出了不同条件下的静力平衡和非静力平衡模式的假设方案 [11]

《2.2 模式垂直坐标与参考大气的选择》

2.2 模式垂直坐标与参考大气的选择

1) 垂直坐标。 由于地球表面的不规则分布，自然高度坐标下的模式下边界很复杂，难以给定下边界条件。GRAPES 模式采用“ 高度地形追随坐标”[12]：。 因为该坐标既有贴地“ 光滑” 坐标特性，有利于模式贴地面层的输送通量计算，同时，该坐标又兼有高度的性质，坐标面不随时间变化，使坐标面上描述的运动在时间上是前后一致的。

2) 模式大气参考廓线。 对垂直运动方程（3）的量纲分析可知，方程右边的第一项（ 气压梯度力项）和第二项（重力加速度项）的量级为大项，引入“模式参考大气：和( 和示参考大气廓线；表示偏离参考大气状态的扰动量），消除垂直运动方程中满足静力平衡的分量，使垂直运动方程的差分计算中重力加速度与气压梯度力之间由“大量级项”变为“ 扰动小量级项”，使之降低与方程中其他项的“量级差”，从而有效地提高垂直运动方程的计算精度。

《2.3 时间离散化方案》

2.3 时间离散化方案

GRAPES 模式采用半隐式－半拉格朗日时间差分方案，该方案已被世界各大业务气象中心的业务数值天气预报模式所广泛采用，如 ECMWF、法国、英国、德国、加拿大、日本、澳大利亚等国家气象业务中心的业务数值预报模式，其计算公式可简化为：

式( 9 ) 中， ψ代表 等模式预报变 量；下标“ D ”表示沿拉格朗日轨迹上游出发点的变量值；右边表示两时间层的线性项 Lψ 和非线性项 Nψ 权重平均形式，要求： = 1 。较之欧拉显式时间差分方案或者时间分裂差分方案，半隐式－半拉格朗日时间差分方案[13]是绝对稳定的，理论上可以取很长的时间步长，可大大提高时间差分方案的计算效率，而且具有良好的频散特性[14] 。

对于标量（ θ 和 ），可以将热力学方程和连续方程直接写成式（9）的时间离散形式。 对于矢量场（u，v 和 w） 的时间离散化，需要采用不同于标量场的时间离散化处理。 矢量离散法是目前被认为较好的方法[15]。 采用非线性内插三维质点位移速度的方法计算拉格朗日质点运动轨迹，确定上游出发点。

《2.4 空间离散化方案》

2.4 空间离散化方案

1) 水平网格。 作为一个区域与全球统一的数值模式，经/纬 度格点模式是较理想的 选择。 一方面，经/纬度网格模式的网格设计简单，另一方面，可以很方便地设计有限区域与全球范围的嵌套模式（如英国气象局的统一模式），也容易与海洋格点模式相耦合。 因此，GRAPES 模式水平方向采用等经－纬度的水平网格（ 见图 1） 和 格式[16]变量跳点分布设置（见图 2），以及二阶精度的中央差分格式。

《图1》

图1 全球与区域统一模式的经纬度 网格系统设计[3，4]

Fig.1 Design of latitude-longitude grid[3，4] for the global-regional unified model

《图2》

图2 水平交叉 Arakawa-C 网格分布[3，4]

Fig.2 Arakawa-C staggered grid[3，4]

2) 垂直离散化｡ 静力平衡模式由于忽略了垂直运动预报方程( 3 ) ，多采用 Lorenz 的跳层方法设置模式预报变量的垂直分布（见图 3）[17] 。而对于采用三维运动方程的非静力平衡模式来说，采用非均匀的 Charney-Philips 跳层方法设置模式预报变量的垂直分布[18]，有利于减少从模式动力框架计算到物理过程参数化计算的温度场垂直插值计算，减小了垂直气压梯度力项的计算误差，并且可以消除虚假的斜压不稳定[19] 。有研究指出[20]，对于处理强涡流与强层流系统、地转适应过程等问题，CharneyPhilips 的跳层[18] 设置优于 Lorenz 的跳层[17] 设置。目前，英国气象局的新一代模式采用 Charney-Philips 的跳层设置[18] 。因此，GRAPES 模式采用 CharneyPhilips 的跳层[18] 设置。

《图3》

图3 变量的 Charney-Phillips 和 Lorenz 跳层分布

Fig.3 Vertical arrangement of the model prognostic variables with Charney-Phillips’ and Lorenz’s methods

注：分别表示模式预报变量垂直速度、位温、Exner 气压函数和东西、南北方向的水平风

《2.5 模式程序软件的设计策略》

2.5 模式程序软件的设计策略

为了搭建一个多尺度通用的或者全球／区域一体化的数值天气预报模式，除了网格系统的选择、模式大气的开关设置（δ1､δ2）以外，还要从软件工程构造上采取适当的策略。参考国外其他国家数值预报模式程序软件的设计策略，尤其是天气预报模式（WRF）的程序软件包基础架构，设计构建一套全新的、统一的 GRAPES 模式程序软件包，达到一套模式程序软件包，满足多种研究与业务预报的应用需求。

1) 标准化｡参照其他领域的软件标准和美国、ECMWF、法国、英国等的数值天气预报系统软件编程标准，建立了一套既要与国际接轨又要满足开发合作与持续创新的实际需要的、基于FORTRAN-90 的气象数值预报模式程序编写标准。另外，建立一套辅助程序转换实用工具，将已有的 FORTRAN-77 程序半自动转换为 FORTRAN-90 程序，以弥补采用手工操作的程序移植方式引起的易出错、修改不全面等的不足。

2) 模块化｡ 所谓模式程序设计的模块化（在 FORTRAN-90 计算机语言中，称为“Module”）就是将数值天气预报系统中的科学计算内容，按照不同的数值计算方法、计算内容、计算顺序、不同的物理过程参数化方案等分为独立计算的模块进行编程，达到根据不同的需要“插拔式”地选取不同的模块，组装成针对不同的预报对象、不同研究目的的数值预报模式，达到一套模式程序软件包却有多种研究与业务预报用途的目的。例如，全球中期数值天气预报模式、有限区域短期数值天气预报模式、中尺度数值天气预报模式、热带台风数值预报模式均是由同一套“模块化”的程序软件包组装而成。

3) 并行化｡多节点多处理器体系结构的并行计算机是未来高性能计算机的主流发展方向。因此，要使数值天气预报系统在高性能并行计算机环境下高效运行，对其模式程序作并行化计算处理是不可缺少的重要工作。现有并行化数值天气预报模式程序存在两方面的问题，一是大部分气象科学家编写的数值模式，其并行化计算效率不高；二是经并行化计算编程专家优化的数值模式，其并行化计算效率高，但模式程序的可读性较差。而且，由于高性能并行计算机的快速更换，数值模式中并行化计算的程序需要作相应的大量修改、优化。为了更好地发挥气象科学家和并行化计算编程专家的优势作用，并同时考虑到高性能并行计算机快速更换所带来的影响，笔者研究团队设计了将模式程序分为“驱动层”、“中介层”和“模式层”的“三层结构”的模式程序软件体系（见图 4），使并行优化与串行模式编程既可独立进行，又可有机组合，还可使程序易读、易写、易改、易维护，可移植性强，通用性和灵活性高。如图 4 所示，顶层为驱动层，该层与并行计算机相关，底层为模式层，该层程序基本上是串行程序，是大部分气象科学家的工作层。在顶层和底层之间为中介层，将驱动层与模式层有机地连接起来，在模式程序编译和运行时实现全套模式在负载较平行的条件下并行化高效计算。

《图4》

图4 GRAPES 模式程序软件分层结构

Fig.4 The layer-structure of the GRAPES model code software

《3 GRAPES 模式的应用现状》

3 GRAPES 模式的应用现状

自 2001 年以来,在国家“十五”科技攻关计划､“十一五”科技支撑计划的支持下, 中国气象局持续投入, 坚持自主研究发展了新一代数值天气预报业务系统 GRAPES ｡ GRAPES 模式有两个系列的应用版本 : GRAPES_Meso ( 区域中尺度模式 ) 和GRAPES_GFS ( 全球中期预报模式 ) ｡ 各区域气象业务中心又在 GRAPES_Meso 模式的基础上, 发展针对本区域天气特点的业务与研究模式, 如广州热带海洋气象研究所 ( 简称广州热带所 ) 的 GRAPES_TMM ( 热带气象模式 ), 上海台风研究所 ( 简称上海台风所 ) 的 GRAPES_TCM ( 热带气旋模式 ), 成都高原气象研究所 ( 简称成都高原所 ) 的 GRAPES_TPM ( 高原气象模式 ) 等｡ 图 5 列举了 GRAPES_Meso 和 GRAPES_GFS 的主要发展历程｡ 2001 年 1 月 , GRAPES 计划正式启动｡ 2004 年 3月, GRAPES_Meso V1.0 ( 60 km ) 第一版建立, 并在国家气象中心､ 广州热带所实时试运行应用｡ 2005 年 5 月, GRAPES_Meso 升级为 V1.5 ( 30 km ), 在国家气象中心､ 广州热带所､ 上海台风所实时运行应用｡ 2006 年 7 月, GRAPES_GFS ( 110 km ) 全球中期预报系统建立｡ 2008—2010 年, GRAPES_Meso 继续改进升级为 V2.5 ( 150 km ) 版本, 在国家气象中心､广州热带所､上海台风所､成都高原所业务运行应用｡ GRAPES_GFS ( 110 km ) 全球中期预报系统正式版本确认 V1.0 ( 55 km ), 并在国家气象中心业务运行应用｡ 2011 年至今, GRAPES 区域模式和全球模式继续在改进升级､应用｡ 表 3 给出了广州热带所的 GRAPES_TMM 模式对南海热带气旋 ( 台风 ) 路径预报误差的结果｡ 从表 3 可以看出, 路径预报误差逐年在减小, 也即预报水平在提高, 对气象减灾防灾保障服务起到了很好的科技支撑作用｡

《图5》

图5 GRAPES 模式发展历程

Fig.5 The historic millstones of GRAPES development

《表3》

表3 南海热带气旋（台风）路径预报误差

Table 3 The distance error of tropical cyclones over South China Sea by GRAPES model

注：数据由广州热带所提供

GRAPES 模式除了在气象部门的业务中心得到应用以外，在大气科学研究与应用培训工作、集合预报、大气科学观测模拟试验、沙尘暴预报、水文洪水预警预报、热带气旋预报、快速更新循环（RUC）等领域也应用十分广泛，未来还将向气候模拟、闪电预报等方向延伸（见图 6）。

《图6》

图6 GRAPES 模式的主要应用领域

Fig.6 The main applications of GRAPES

注：实线框表示已形成业务能力，虚线框表示尚未形成业务能力，处于研究阶段

《4 GRAPES 模式的发展》

4 GRAPES 模式的发展

GRAPES 模式的进一步发展涉及模式初值化、模式动力框架、模式物理过程、专业应用模式拓展、模式不确定性与集合预报等。 这里只就模式动力框架中的垂直坐标和网格系统问题进行讨论。

《4.1 模式垂直坐标》

4.1 模式垂直坐标

为了易于处理模式计算域的下边界条件，气象数值预报模式多采用下边界光滑的地形追随坐标， 其模式面随地形高低起伏，使得水平气压梯度力的计算由原来的一项差分计算变成了两项差分计算之差，其中一项是与地形坡度相联系的项，它会导致水平气压梯度力计算的误差，这种误差影响从模式低层一直到模式顶层。而且，随着模式水平分辨率的提高，模式地形坡度越“陡峭”，引起的这类计算误差越大。对此，有几种处理方法。 第一种是修改气压梯度力项的计算方案，例如：回插等压面法、递推算法[21]、扣除法[22,23]、基于静力方程订正的气压回插法[24]等｡第二种是 η- 阶梯地形坐标[25]，即用等高面将模式地形切割成“台阶”状，水平等高面就是模式面。第三种是构造新的地形追随坐标面随高度衰减的函数，加快“弯曲”的地形追随坐标面随模式高度增大而逐渐变水平[26] 。

第三种方法的优点是通过对衰减系数的调节，使得垂直坐标在底层受到地形作用的影响逐渐减小，弯曲的坐标面逐渐趋缓，到了高层坐标面趋于平直，更接近大气的真实状况（见图 7）。这一点可以从气压梯度力转换计算公式加以理解。

《图7》

图7 不同高度的水平模式面

Fig.7 Vertical level of model surface

式（10） 左边项为自然高度坐标下的气压梯度力项，经地形追随坐标转换后变成等式右边的两项，其中右边第二项是和地形坡度（b） 以及地形坡度雅克比（ ）有关的项， 随高度减小，至高层趋于零，这时等式右边第一项等于左边项，意味着高层大气运动趋于水平。 GRAPES 模式将参考第三种方法对垂直坐标加以改进。

《4.2 模式水平网格系统》

4.2 模式水平网格系统

经纬度网格是世界上大多数国家气象业务中心全球业务数值预报模式最常用的水平网格系统，其优点在于容易构造应用，而且其正交的经纬线可以和南北 －东西风方向相一致，更易于描述模式大气的水平运动和平流输送计算。 但是，经纬度水平网格系统存在致命的缺点。 经线 在南－北两极辐合（见图 1），使得靠近极区的格点距离要大大小于中低赤道地区的格点距离（ 如模式水平分辨率取 0.1 °×0.1 °时，在赤道附近的格距约为 11.1 km，而在极点附近的格距仅为 0.02 km， 两者 之比达 573） ，并最终在两极出现奇异点。 这些给模式的差分计算带来了极大的麻烦，解决的办法如下。

1) 精简格点法。 即将靠近极区处沿纬圈的格点数减少，而且越靠近极点的纬圈，沿纬圈的格点数精简越多。 大多数谱展开模式采用此方法。

2) 球冠投影法。 即将两极的区域（“ 球冠”）投影到一个均匀网格面上进行差分计算，然后再映射回两极区域。 大多数格点差分模式采用此方法。

3) 设计一个分布较均匀的水平网格系统，避开传统的经纬度网格经线辐合引起的计算不稳定的问题，如三角形网格、六角形网格、多边形网格等。 此外，一种灵感可能来自中国古老的道教阴阳八卦理念的新网格系统已被提了出来，这就是所谓的“阴阳” 网格系统[27]，或称叠交网格系统（见图 8c）。

《图8》

图8 阴阳网格设计

Fig.8 Design of a Yin-Yang grid

阴阳网格系统由两个相似的经纬度有限区域网格构造成一个完整的全球网格系统，该网格具有格点分布准均匀和“对称”、无奇异点的优点。 该网格的提出为有效解决经纬度网格在近极区格点辐合及极地奇异点计算问题提 供了新的思路和途径。 另外，阴阳网格系统的基础仍然是经纬度网格，仅仅是一个坐标旋转变换的问题。

式（11）中，， ，上标“ i” 表示阴网 格；上标 “ａ” 表示阳网格；球坐标的计算域为：

式（12）中，δ 为设定的阴阳网格重叠区。因此，很容易地将经纬度系统转换成阴阳网格系统，自然也很容易继承原来经纬度网格的程序软件｡ Peng 等于 2006 年成功地实现了在阴阳网格系统上的欧拉平流方案的理论计算[28] 。Li 等也于同年成功地实现了在阴阳网格系统上的拉格朗日平流计算[29] 。GRAPES 模式将采用阴阳网格系统以解决传统经纬度网格的极区计算问题，并尝试全球－区域的双重（或多重）同步双向嵌套的一体化运行策略。

《5 结语》

5 结语

我国科学家自主研究发展的 GRAPES 模式在科学设计上采取了单一的动力框架、多种物理过程参数化方案可选的科学策略，在软件工程上采用模块化、标准化、并行化的软件程序编码技术，使之成为多尺度统一模式，可以用于中尺度、短期区域、全球中期天气预报，区域与全球气候预测，沙尘暴气溶胶环境模拟等。水平尺度可从 1 km 到 100 km，时间尺度可从几小时到 2 个星期，甚至到月、季、年。这一模式科学设计所带来的优点是显然易见的。

1) 使科研开发与业务预报模式平台统一，可加快科研成果向业务应用的转化。

2) 可避免一个业务中心同时运行多套不同构架的模式系统的局面，大大降低了业务运行及维护成本，包括数值预报模式系统在不同体系结构的计算机上的移植成本。

3) 便于更多领域（如数值模式动力框架、资料 同化、物理过程参数化、软件工程等） 的研发人员参与，集约化发展。

作为一个完整的数值天气预报系统，模式积分的初始值生成方案———资料同化系统也是至关重要的。GRAPES 的资料同化系统采用了基于泛函理论的三维／四维变分同化技术（3D/4DVAR-DA）[30]，以便更有利于同化应用更多雷达、卫星遥感观测资料，构造更协调的模式积分初值。变分资料同化技术是目前世界上先进的业务应用技术。

天气气候无缝隙预报将是未来气象数值预报的发展方向。针对天气气候无缝隙预报的问题，Shukla[31] 指出没有科学根据可以人为地画出一条边界来，区分中尺度预报、天气预报、季节预测、ENSO（厄尔尼诺和南方涛动的合称）预测，年代际预测以及气候变化。但是，从计算机资源、模式复杂程度的实际情况考虑，可以建立满足不同时间尺度需要的天气气候预报预测系统，而从研究的角度则可以在一个统一的框架下发展天气气候预报预测系统。GRAPES 模式是多尺度通用的一体化模式，也为未来我国的天气气候无缝隙预报预测系统的研究发展奠定了基础。