下一代大气模式中的数值方法综述

沈学顺，李兴良，陈春刚，唐杰

(1.中国气象局地球系统数值预报中心，北京 100081；
2.中国气象科学研究院灾害天气国家重点实验室，北京 100081；
3.西安交通大学航天航空学院，陕西 西安 710049；
4.西安交通大学航天航空学院机械结构强度与振动国家重点实验室，陕西 西安 710049)

数值天气预报的萌芽可以追溯到1904年挪威科学家Vilhelm Bjerknes提出把天气预报看成求解描述大气运动控制方程组的科学思想，即通过实时观测的气象数据求解控制方程组来预测未来时刻的天气状态。实践证明，解析求解这些控制方程组的天气预报方法是行不通的。1922年，英国气象学家RICHARDSON[1]创造性地用数值方法求解大气运动方程组，开创了用数值方法预报天气的革命性技术，即数值天气预报。随着数值算法、高性能计算机技术、气象探测技术以及动力气象学、数值预报理论的不断发展，数值天气预报在近100年时间内取得了巨大的成功[2-3]。数值天气预报已成为现代天气预报业务的核心支撑，在气象防灾减灾、社会公众服务、国民经济建设和国防等工作中发挥着不可替代的作用。

数值天气预报系统由大气数值模式、资料同化及模式后处理等部分构成，其中数值模式是核心。大气数值模式主要由动力框架和物理参数化过程两大模块组成。大气模式中最早采用的数值方法是基于规则网格点的有限差分方法，该方法通过差分算子近似大气控制方程组中的导数项[4-6]。随后，自20世纪70年代末开始，谱模式在大气模式中占主导地位，并在业务模式中应用最为广泛[7]。谱方法(spectral method)通过一系列正交基函数对控制方程中预报量进行逼近，实现空间离散。它相比于早期格点模式的优势在于能有效避免球面的极区问题和消除相速度误差(phase-speed error)等[8]。随着模式向高分辨率方向发展，谱模式勒让德变换在高分辨率下的效率瓶颈和难以在众核超级计算机上实现高效并行化计算的缺点日益突出，基于格点模式的数值方法(如：半拉格朗日方法(semi-Lagrangain method)[9]和有限体积法(finite volume method)[10]等)在进入21世纪后再次兴起。

近年来，随着异构众核高性能并行计算机和E级计算的高速发展，大气数值模式中涌现出一系列能够高效利用超级计算机性能和保证良好计算效率的高扩展算法，高扩展算法已成为下一代大气数值模式的研究焦点[11]。世界上各主要业务中心已于2010年左右开始下一代高可扩展模式研发，如：欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts，ECMWF)于2015年开始的用于E级计算的天气预报高能效规模算法(Energy-efficient Scalable Algorithms for Weather Prediction at Exascale，ESCAPE)计划，英国气象局于2011年开始的“工合(GungHo)”(Globally Uniform Next Generation Highly Optimized)计划，以及美国国家海洋和大气管理局(National Oceanic and Atmospheric Administration，NOAA)于2016年开始的下一代全球预报系统(Next Generation Global Prediction System，NGGPS)计划等。

一般而言，要实现未来大气状态的气象要素预报，大气数值模式需要基于球面计算网格对大气动力与热力学控制方程进行空间数值离散，继而采用时间离散格式向前推进模式积分。为此，本文将从数值算法、准均匀球面网格和时间积分方法等3个方面对国内外可扩展的下一代大气模式进行综述，希望对相关领域同行的研究起到帮助和推动作用。需要指出的是，全球模式与区域模式的差异体现在模式是否有边界处理，而在数值算法、球面网格和积分方案方面并无显著区别，因而本文所介绍内容对全球模式与区域模式均适用。

自20世纪20年代RICHARDSON[1]的创始性工作以来，应用在大气模式中的数值方法呈现多样性，如传统格点模式采用的有限差分方法、有限体积方法、半拉格朗日方法[12-13]和非格点模式(ECWMF)采用的谱方法等。2007年，WILLIAMSON[7]在其综述“全球大气模式动力框架演变”中详细回顾了传统数值方法在全球大气模式中的应用进展。

近年来，高性能并行计算机的迅速发展对应用于大气模拟的数值计算方法提出了新的挑战，即：数值计算方法必须能够充分利用超级计算机所提供的计算潜力[14-15]。基于这一背景，过去十多年中多种适用于大气数值模拟并能高效利用超级计算机的高精度数值方法被应用于发展新一代大气模式，代表性的方法包括：谱元(spectral element)方法[16-20]、间断伽辽金(discontinuous Galerkin)方法[21-25]、混合有限元(mixed finite element)方法[26-29]以及多矩约束有限体积(multi-moment constrained finite volume)方法[30-39]。作为有潜力的下一代大气数值模式核心基础算法，这些数值方法具有高精度、守恒、网格适应性强和可扩展等特点，国际先进业务研究中心开始或已经采用它们来发展新一代大气数值模式。与传统的有限差分和有限体积数值方法不同，这些数值方法能够在局地单元开展高阶多项式数值重构，因而获得了高精度和良好可扩展特性。然而，上述高精度格式其计算过程较为复杂、数值计算鲁棒性不高、需要与合适的时间积分方案相配合等，这些要求给大气模式开发者采用这类高精度数值方法设计下一代数值模式提出了新的挑战。

1.1 有限元方法

有限元方法最早发展于20世纪40年代初[40]，并在50年代后期取得重要进展并得到广泛应用。有限元方法的空间离散格式一般通过加权余量法(weighted residual method)导出。在实际应用中，通常采用伽辽金(Galerkin)方法，即：选取离散空间上的基函数(有时又称为试探函数(test function))作为权重函数，使之与误差函数具有正交性，得到空间离散公式。它的优势在于能够灵活地处理复杂几何剖分。

有限元方法在固体力学领域得到了广泛应用[41]，在流体计算中主要应用于低马赫数流动的模拟[42]。多种有限元方法被应用于地球流体力学数值模拟，如：基于有限元方法的海洋模式[43-48]。有限元方法应用到大气模拟可追溯到20世纪60年代[49]，这一领域至20世纪80年代涌现出更多工作。STANIFORTH[50]对有限元方法在大气模拟中的应用做过综述，强调有限元模式有望更好地探索、理解和预报潜在更复杂的天气现象。BÉLAND et al.[51]采用有限元方法研究了罗斯贝(Rossby)和重力正规模态(normal modes)模拟的准确性，并且与二阶有限差分垂直离散方案做了比较。

相比较于计算域全局展开的谱方法，有限元方法仅在局部单元进行有限元函数级数展开，没有全局并行通信，因而具有良好的计算效率。然而，与有限差分方法相比，有限元方法空间数值离散后需要求解代数方程组，因此它的复杂程度更高，计算代价相对较大[52]。此外，由于单元内未知量多于一个，有限元方法存在潜在计算模(computational mode)[53]问题，该计算模是非物理解。特别地，当表征小尺度大气波动时，非线性项的数值离散容易引起错误的数值解。为了抑制计算模，有限元方法通常需要选择性地采用耗散机制[54]，这可能是该方法在大气模式水平离散中没有得到广泛应用的重要原因。相反，采用有限元方法作为垂直离散格式在大气模式中得到良好应用，因为大气数值模式的垂直数值耗散对有限元方法抑制计算模发挥了关键作用。BURRIDGE et al.[55]在欧洲中期天气预报中心(ECMWF)的预报模式中使用有限元方法完成了垂直离散；
UNTCH and HORTAL[56]在ECMWF模式中引入垂直有限元平流离散方案，相对于有限差分模式极大提高了线性重力波位相速度模拟精度，减小了洛伦兹(Lorenz)交叉跳点的计算模和垂直计算噪音。

英国气象局下一代数值天气预报模式动力框架“工合(GungHo)”研究计划将混合有限元方法作为数值方法的首选。STANIFORTH et al.[28]研究表明，混合有限元离散中速度自由度与气压自由度数量比值小于2时，会出现明显的虚假惯性重力波，而两者比值大于2时，则发现类似于六边形网格上的虚假罗斯贝波[57]。COTTER and SHIPTON[26]针对线性浅水波方程进行了研究，表明混合有限元方法在三角形网格、六边形网格和正方形网格上通过选择适合的基函数对(element pair)，调整气压和速度自由度数量比率，进而可以避免虚假的地转模[58]，如：重力惯性模和罗斯贝模。此研究发现，在三角网格中对速度离散采用一阶Brezzi-Douglas-Fortin-Marini空间而气压采用线性不连续空间，在四边形网格中速度离散选用k阶Raviart-Thomas空间而气压则用k阶不连续空间，就可以保证速度自由度和气压自由度数量之比为2∶1，从而避免虚假的地转模，并最小化频散和耗散误差。同时，对非线性系统连续方程采用这样的有限元空间并运用间断伽辽金方法可以保证质量守恒。

实际上，混合有限元方法[59]是跳点网格(如C网格)上有限差分方法的拓展，对于不同变量选用不同有限元空间函数(基函数)，规避了有限差分方法在A格点上的虚假气压模，同时不受数值网格正交性的约束。目前，这类应用混合有限元方法在球面浅水波模拟中已获得良好效果[28-29,60-61]，进一步推广至三维完全可压缩大气模式及针对不同变量的基函数选择还在进一步研究当中[59]。当前，MELVIN et al.[62]基于混合有限元离散格式在笛卡尔坐标下发展了三维非静力动力框架，它不仅在模式精度上与现有英国气象局业务模式ENDGame动力框架相当，而且具有潜在的众核高性能计算与不受模式数值网格约束的优势。

1.2 谱元方法

MA[63]最早将谱元方法应用于地球流体力学模拟，建立了适用于复杂几何域的谱元浅水模式，研究表明由于谱元方法的耗散和频散误差小，模式能够准确地模拟海洋锋面现象。TAYLOR et al.[64]将谱元方法应用于球面浅水模式，评估了谱元方法用于气候模式的可行性，讨论了谱元模式在气候模拟中的潜在优缺点。FOURNIER et al.[19]将谱元方法应用到静力原始方程，建立谱元大气模式(Spectral Element Atmospheric Model，SEAM)。随着当时大规模分布式高性能计算机的出现，谱元方法的高可扩展性得到了印证。GIRALDO and ROSMOND[18]在水平方向采用谱元方法发展了一套新的数值天气预报动力框架，测试了三个正压三维标准测试和四个斜压三维标准测试，如：四波罗斯贝波等；
并与四个天气预报模式和气候模式相比较，模拟结果表明谱元大气模式和谱变换模式精度相当，该模式在分布式计算机上的并行度可实现线性可扩展。THOMAS and LOFT[65]建立了基于立方球网格的动力框架，该框架在水平方向使用高阶谱元方法，采用垂直混合气压坐标和半隐式时间积分方案。

当前，DENNIS et al.[15]发展的CAM-SE模式、美国海军研究院具备谱元可调的非静力通用大气模式(Nonhydrostatic Unified Model for the Atmosphere，NUMA)以及公共地球系统模式(Community Earth System Model，CESM)等均是采用谱元方法发展的新一代大气模式。韩国大气预报系统研究院(Korea Institute of Atmospheric Prediction Systems，KIAPS)在韩国气象厅下一代大气模式中采用高精度谱元方法，具体设计方案是在水平方向采用谱元方法离散，在垂直方向采用有限差分离散方法[66-67]。谱元方法具有良好的局地计算密集型特征，在分布式高性能计算机上可获得高可扩展性[15-16]，同时具有谱收敛和几何灵活(geometrical flexibility)等优势，已成为下一代预报模式的优选数值计算方法之一。

1.3 间断伽辽金方法

间断伽辽金方法，与前述有限元和谱元方法不同，数值解和基函数均局限于单元网格内，相邻网格内的近似函数和权重函数在单元界面上的值一般不相等。对于守恒方程，相邻单元网格界面的数值通量一般通过求解(近似)黎曼(Riemann)问题得到，因此有限体积方法中发展的(近似)黎曼问题求解器均可用于间断伽辽金方法。考虑这一特征，间断伽辽金方法也被认为是一种有限元和有限体积相混合的方法。间断伽辽金方法和谱元方法一样也采用不等间距分布的节点，但其数值解定义点不包括单元边界，如采用Legendre-Gauss(LG)点可以获得非常高的计算精度。由于间断伽辽金方法中引入了有限体积的概念，它可以保证数值守恒性。

REED and HILL[68]最早引入间断伽辽金方法用于求解静态中子线性平流问题。LESAINT and RAVIART[69]对间断伽辽金方法应用于线性平流方程的数值特性进行了分析。CHAVENT and SALZANO[70]首次将间断伽辽金方法应用到求解非线性方程，研究发现若将一阶多项式间断伽辽金方法空间离散和简单向前时间积分相结合一般计算不稳定，仅在特定情况下可保证计算稳定性。COCKBURN and SHU[71-72]采用TVB(total variation bounded)限制器、龙格-库塔(Runge-Kutta，RK)时间积分和间断伽辽金方法求解非线性方程获得了极大成功。随后涌现了大量采用间断伽辽金方法来求解实际问题的研究工作。SCHWANENBERG and KÖNGETER[73]采用龙格-库塔间断伽辽金方法研究了一维浅水方程求解及其在水利上的应用，充分展现了其质量和动量守恒、几何灵活和紧致性(compactness)好等优点。COCKBURN and SHU[74]对采用龙格-库塔时间积分和间断伽辽金方法发展的数值模型做了较详细的综述。

在大气模式领域，GIRALDO et al.[75]基于正二十面体网格在笛卡尔坐标下发展了间断伽辽金浅水波模式，其中空间离散采用nodal类型展开。NAIR et al.[76-77]基于立方球网格采用modal类型展开的方式发展了间断伽辽金平流模式和浅水波模式，研究结果表明间断伽辽金方法的精度与标准谱转换模式相当。GIRALDO and RESTELLI[22]将间断伽辽金方法应用于求解通量形式的Navier-Stokes方程，研究了非静力中尺度大气现象，数值试验表明：在热泡和山脉波试验中间断伽辽金非静力模式表现更好；
考虑计算效率，在最坏情况下，间断伽辽金模式的运行速度要比非守恒谱元模式慢近一半，在最好情况下，间断伽辽金模式的运行速度与守恒的谱元模式相当。尽管目前采用间断伽辽金方法建立的大气模式还未用于实际业务预报，然而基于间断伽辽金方法的研究型大气模式一直在发展中，比如美国国家大气研究中心(National Center for Atmospheric Research，NCAR)的HOMME(high-order method modeling environment)静力模式[24]，美国海军研究院的间断伽辽金全球、区域一体化NUMA[78]，德国的DUNE(Distributed and Unified Numerics Environment)大气模式[79]等。

1.4 多矩约束有限体积方法

有限体积方法是基于物理守恒律发展的一种数值方法，它能够严格保证数值解的守恒属性。近年来人们对具有良好局地守恒属性，如紧致模板、高强度浮点计算和极少数据通信的数值方法产生新的兴趣。为此，设计能够适应众核超级计算机的一类有限体积方法在非静力大气模式[35,80-83]中的应用重新受到科学家们的关注和研究。

作为对传统有限体积法的扩展，多矩约束有限体积法是一种数值求解流体力学方程组的新方法。多矩方法借鉴了CIP(constrained interpolation profile)方法的基本思想[84-87]，在算法构造中同时采用不同类型的矩，如：积分平均值、状态变量的点值及其导数值，不同矩可以采用不同形式的控制方程进行预报，这些预报方程必须保持与原守恒律方程相容。它的构建过程基于明确的物理概念，具有良好可扩展性、可达到高精度、并能保证严格的数值守恒。与间断伽辽金和谱元方法相比较，它较好地平衡了数值模式在计算精度、鲁棒性、算法简洁性和计算效率等各方面的要求，具有较高的实用价值。不同于传统有限体积方法在单元网格内定义单个自由度，多矩约束有限体积方法与谱元和间断伽辽金方法类似，在单元网格内定义多个自由度(或者未知量)，且自由度定义节点可以自由选择，包括等距节点和Legendre-Gauss-Lobatto(LGL)节点等。该方法通过施加不同类型的约束条件，如：物理量的积分平均值、点值及导数值，导出自由度的演变方程，最后再通过龙格-库塔时间积分获得未知量(或自由度)在新时刻的预报值。由于多矩方法在空间拉格朗日(Lagrange)插值重构中引入了多矩约束条件，相比有限元方法中的等距拉格朗日插值具备更好的计算稳定性，能有效抑制高阶重构函数在网格单元边界处的数值振荡[32]。XIAO et al.[34]发展了一套基于多矩约束通量重构的通用高精度数值架构，它可以解释成拉格朗日和埃尔米特(Hermite)混合插值方法，可以从更一般的途径构造高精度数值格式。同时研究[88]表明，多矩约束方法在相同自由度下计算稳定的最大库朗条件数(Courant-Friedrichs-Lewy(CFL)数)要大于间断伽辽金方法。

多矩约束方法已被成功用于发展地球流体力学数值模式。CHEN and XIAO[33]采用多矩约束有限体积方法在立方球面网格上建立了全球浅水波模式。该模式中状态变量的点值采用微分形式方程组通过求解导数黎曼问题来更新，积分平均值采用通量形式方程组用有限体积公式更新保证数值守恒性。LI et al.[89]采用球面阴阳网格发展了球面浅水波模式，不同的是状态变量点值采用半隐式半拉格朗日方法来更新。尽管对状态变量点值的更新方法有所不同，但这两个浅水模式对球面大气波动，如：罗斯贝波、过山气流波等，均获得良好的数值模拟结果，可与国际其他先进模式相媲美。采用多矩约束有限体积框架，可以构造基于局地重构的任意高阶格式。采用三阶、四阶和五阶多矩约束有限体积方法，基于不同球面准均匀网格，如：阴阳网格、立方球网格、正二十面体三角形网格和六边形网格，已发展了各类全球浅水模式[36,38-39,90-91]。综合考虑模式的计算精度、效率等因素，多矩方法相比其他先进数值方法有其独特优势。

LI et al.[35]首次将多矩约束有限体积方法应用于完全可压缩、非静力大气模式。一系列标准数值试验，如：密度流试验、惯性重力波试验、热泡试验、Schar复杂地形试验以及线性静力/非静力钟型山脉波试验[37]等，均表明该模式的计算精度与当前国际先进谱元模式和间断伽辽金模式相当。多矩约束有限体积方法在发展球面浅水模式和非静力大气模式中获得了良好数值结果，有望为发展我国下一代高精度、高可扩展、守恒的天气模式甚至气候模式提供可靠的动力学数值架构。

在大气数值模式中代表性的球面网格有：经纬度网格、阴阳网格、立方球网格和正二十面体网格(图1)。几乎所有的业务中心均采用过经纬度网格，它具备良好的正交性，方便进行差分计算。然而在超级高性能计算环境中，经纬度网格的极区问题及其带来的数据通信瓶颈变成了制约全球模式性能的一个关键因素。

图1 大气模式中代表性球面网格示意图(a.经纬度网格，b.阴阳网格，c.立方球网格，d.正二十面体网格)Fig.1 Schematic interpretation of representative spherical grid in atmospheric model (a. latitude-longitude grid, b. Yin-Yang grid, c. cubed-sphere grid, d. icosahedral grid)

KAGEYAMA and SATO[92]提出了准均匀、无数值奇异的阴阳网格，它由两个低纬度经纬度网格相互扣接而成。阴(或阳)网格是普通经纬度网格一部分，继承了普通经纬度网格特性，且低纬度网格单元面积较均匀，并避免了普通经纬度网格数值计算奇异性问题。球面阴阳网格属于组合网格，在阴阳边界处的数据交换需要通过插值运算实现，为保证质量守恒需要特别的考虑[93]。立方球网格由SADOURNY[94]提出，该网格通过球心投影正立方体表面至球面而成，在其六个投影面上可生成网格单元面积相当均匀的结构化网格。根据投影方式不同，它可分为心射切面投影(gnomonic)立方球网格和保角(conformal)投影立方球网格[95]。前者生成的球面网格不正交，但拥有解析的坐标变换关系；
后者生成的球面网格正交但不存在解析坐标变换关系。SADOURNY et al.[96]提出了正二十面体球面网格，它通过投影正二十面体至球面得到二十个球面三角，而在每个三角内又细分成多个小三角或六边形的网格单元，属于非结构网格。此外，ECWMF从2016年起在业务模式中采用立方八面体网格，它是简化的高斯网格(reduced Gaussian grid)[97]的一种，即八面体简化的高斯网格(octahedral reduced Gaussian grid)，经向格点数(latitude points)保持和简化的高斯网格一样，而相对于简化的高斯网格，其纬向格点数(longitude points)越往极区格点数越少。

网格的几何属性在数值模拟中起到非常重要的作用，它和数值方法结合在一起共同影响数值模拟的准确性。STANIFORTH and THUBURN[98]对全球天气、气候模式采用的水平网格有详细综述，特别关注了球面水平网格的计算模(computational modes)和格点印记(grid imprinting)问题。他们认为网格几何属性会直接影响数值模式动力框架一些基本的理想性能，如质量守恒、准确描述大气地转平衡和调整问题、计算模消失或者被很好地控制住的问题、虚假快速传播的Rossby波问题、位势和气压梯度不应产生非物理的涡度源、最小化格点印记等。格点印记是数值模式底层网格结构潜在误差的一个典型现象。数值试验观察发现，采用准均匀球面网格的模式其底层网格结构可能会以噪音或者系统误差的形式在模拟结果中凸现出来，如浅水波试验中立方球网格内嵌立方体八个顶点对应的球面点附近的格点印记[33,38,76-77,99]和阴阳网格边界处格点印记[89,100-102]等。

在当前以及未来超级高性能计算条件下，传统数值方法，如：有限差分、有限体积方法等，配合经纬度网格在数据通信问题上遇到计算瓶颈，很可能会制约它的继续发展。毫无疑问，局地高精度、守恒、高可扩展性的数值方法，如：前述的有限元方法、谱元方法、间断伽辽金方法以及多矩约束有限体积方法，具有很好的网格灵活性，在准均匀网格上采用此类数值方法发展模式动力框架是当前的一个热点研究方向[103]。

为保证良好的计算效率，许多全球业务天气预报模式采用半隐式半拉格朗日时间积分方案，比如：平流过程采用半拉格朗日方法来计算，在上游点计算轨迹不交叉情况下时间步长不受限制，而对于声波、重力波等快波则采用半隐式方法来计算[104]。半隐式模式通过适当的线性化和空间数值离散最终得到一个大规模、稀疏和对角化的矩阵，所求解问题转换成大型代数方程，最终采用高效求解器通过迭代获得变量的数值解。

实际大气模式网格剖分时，其水平和垂直网格距之比远大于1(十倍甚至于百倍)，为保证计算效率大气模式常采用水平显式-垂直隐式(horizontally explicit-vertically implicit， HEVI)的时间积分方案，即在水平方向的大气波动项采用显式积分，而隐式处理垂直传播波动项。相对于全局半隐式求解方法，它克服了垂直方向的时间步长稳定性限制，同时避免了并行计算中的全局通信。HEVI显式分裂(split-explicit)方法处理大气控制方程中不同速度的波动时，将方程中的快变项再分成多个小时间步进行积分，同时在垂直方向上用Crank-Nicolson等隐式方法求解[105]，如WRF模式[106]为采用该方法的代表性大气模式。

大气模式中的显式-隐式(implicit-explicit，IMEX)多时间层方案按照波动的快慢分别采用隐式和显式方案进行处理。为了获得二阶或更高精度，一般会用到多于两个时间层的信息。DURRAN and BLOSSEY[107]用IMEX方法求解静力原始方程和可压缩Boussinesq方程。然而，通过稳定性分析可以发现，多时间层的IMEX方案会出现计算模，因此需要通过过滤计算模以得到满意的数值解[108]。相比较而言，IMEX龙格-库塔(Runge-Kutta，RK)方案可以基于单时间层内多个小时间步获得良好的稳定性和计算精度(二阶及其以上)，即使处理刚性(stiff)物理问题[109-110]，在数学上表现为算子特征值极大。ULLRICH and JABLONOWSKI[111]采用算子分裂(operator-split)Runge-Kutta-Rosenbrock(RKR)方法积分非静力大气模式，在垂直方向采用线性隐式Rosenbrock方法离散，而水平方向采用显式龙格-库塔方法积分，可以获得三阶精度。算子分裂RKR方法实际上是HEVI方法，因而它的积分稳定性最后受限于水平方向网格矩和波动速度。值得一提的是，ULLRICH and JABLONOWSKI[111]控制方程在垂直方向所有项均统一采用Rosenbrock方法隐式求解，没有区分波动传播特点，其代价是会减慢垂直平流。WELLER et al.[112]采用半隐式或HEVI方法，理论分析了多种龙格-库塔IMEX方法求解线性双曲问题的精度和稳定性，并提出了一种新的HEVI-split方法，指出该方案可提高重力波模拟的精度和稳定性。相对于半隐式方法，WELLER et al.[112]研究认为HEVI-split方法在保持计算精度上更有优势。BAO et al.[113]采用基于Strang类型的算子分裂HEVI方法积分间断伽辽金非静力大气模式，在水平方向采用鲁棒性强保型龙格-库塔时间积分，垂直方向的隐式积分释放了受声波约束的时间步长，获得二阶时间收敛精度。

面向未来并行超级计算机，高阶完全显式时间积分方案在计算精度方面无疑是不错选择，然而其时间步长受限于声波的快速传播，当前仍然需要在计算精度和效率上进行权衡。综合考虑各种因素，在全球模式中采用水平方向显式求解大气运动方程而垂直隐式求解快速波动的HEVI方法，并配合上述具有局地高精度特性的数值算法，在未来超级计算机平台上是一套极具吸引力的地球流体数值模式架构。

世界气象组织2015年发布的《Seamless Prediction of the Earth System: From minutes to months》(WMO-No.1156)标志着下一代业务数值预报的发展方向在世界范围内形成了共识。在地球系统科学框架下，发展天气气候一体化耦合数值预报系统以及相应的业务预报服务体系成为世界各主要数值预报业务中心的研究和应用重点。同时，众核高性能计算机和E级计算技术的发展也促使下一代业务数值预报模式在数值算法、网格结构等方面做出新的发展。美国国家环境预报中心和德国气象局近年来投入业务运行的立方球网格有限体积模式(NCEP_GFS-fv3)和基于正二十面体网格的ICON模式是这方面的代表。无疑，在准均匀网格上采用局地高精度、守恒、高可扩展性的数值方法是一个重点和热点发展方向。本着抛砖引玉的目的，本文对近年来大气模式中的数值方法进行了系统的分析总结，希望对国内同行在下一代数值预报研究中起到帮助和推动作用。