沈?qū)W順，李興良，陳春剛，唐杰

(1.中國(guó)氣象局地球系統(tǒng)數(shù)值預(yù)報(bào)中心，北京 100081；2.中國(guó)氣象科學(xué)研究院災(zāi)害天氣國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081；3.西安交通大學(xué)航天航空學(xué)院，陜西 西安 710049；4.西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049)

引言

數(shù)值天氣預(yù)報(bào)的萌芽可以追溯到1904年挪威科學(xué)家Vilhelm Bjerknes提出把天氣預(yù)報(bào)看成求解描述大氣運(yùn)動(dòng)控制方程組的科學(xué)思想，即通過(guò)實(shí)時(shí)觀(guān)測(cè)的氣象數(shù)據(jù)求解控制方程組來(lái)預(yù)測(cè)未來(lái)時(shí)刻的天氣狀態(tài)。實(shí)踐證明，解析求解這些控制方程組的天氣預(yù)報(bào)方法是行不通的。1922年，英國(guó)氣象學(xué)家RICHARDSON[1]創(chuàng)造性地用數(shù)值方法求解大氣運(yùn)動(dòng)方程組，開(kāi)創(chuàng)了用數(shù)值方法預(yù)報(bào)天氣的革命性技術(shù)，即數(shù)值天氣預(yù)報(bào)。隨著數(shù)值算法、高性能計(jì)算機(jī)技術(shù)、氣象探測(cè)技術(shù)以及動(dòng)力氣象學(xué)、數(shù)值預(yù)報(bào)理論的不斷發(fā)展，數(shù)值天氣預(yù)報(bào)在近100年時(shí)間內(nèi)取得了巨大的成功[2-3]。數(shù)值天氣預(yù)報(bào)已成為現(xiàn)代天氣預(yù)報(bào)業(yè)務(wù)的核心支撐，在氣象防災(zāi)減災(zāi)、社會(huì)公眾服務(wù)、國(guó)民經(jīng)濟(jì)建設(shè)和國(guó)防等工作中發(fā)揮著不可替代的作用。

數(shù)值天氣預(yù)報(bào)系統(tǒng)由大氣數(shù)值模式、資料同化及模式后處理等部分構(gòu)成，其中數(shù)值模式是核心。大氣數(shù)值模式主要由動(dòng)力框架和物理參數(shù)化過(guò)程兩大模塊組成。大氣模式中最早采用的數(shù)值方法是基于規(guī)則網(wǎng)格點(diǎn)的有限差分方法，該方法通過(guò)差分算子近似大氣控制方程組中的導(dǎo)數(shù)項(xiàng)[4-6]。隨后，自20世紀(jì)70年代末開(kāi)始，譜模式在大氣模式中占主導(dǎo)地位，并在業(yè)務(wù)模式中應(yīng)用最為廣泛[7]。譜方法(spectral method)通過(guò)一系列正交基函數(shù)對(duì)控制方程中預(yù)報(bào)量進(jìn)行逼近，實(shí)現(xiàn)空間離散。它相比于早期格點(diǎn)模式的優(yōu)勢(shì)在于能有效避免球面的極區(qū)問(wèn)題和消除相速度誤差(phase-speed error)等[8]。隨著模式向高分辨率方向發(fā)展，譜模式勒讓德變換在高分辨率下的效率瓶頸和難以在眾核超級(jí)計(jì)算機(jī)上實(shí)現(xiàn)高效并行化計(jì)算的缺點(diǎn)日益突出，基于格點(diǎn)模式的數(shù)值方法(如：半拉格朗日方法(semi-Lagrangain method)[9]和有限體積法(finite volume method)[10]等)在進(jìn)入21世紀(jì)后再次興起。

近年來(lái)，隨著異構(gòu)眾核高性能并行計(jì)算機(jī)和E級(jí)計(jì)算的高速發(fā)展，大氣數(shù)值模式中涌現(xiàn)出一系列能夠高效利用超級(jí)計(jì)算機(jī)性能和保證良好計(jì)算效率的高擴(kuò)展算法，高擴(kuò)展算法已成為下一代大氣數(shù)值模式的研究焦點(diǎn)[11]。世界上各主要業(yè)務(wù)中心已于2010年左右開(kāi)始下一代高可擴(kuò)展模式研發(fā)，如：歐洲中期天氣預(yù)報(bào)中心(European Centre for Medium-Range Weather Forecasts，ECMWF)于2015年開(kāi)始的用于E級(jí)計(jì)算的天氣預(yù)報(bào)高能效規(guī)模算法(Energy-efficient Scalable Algorithms for Weather Prediction at Exascale，ESCAPE)計(jì)劃，英國(guó)氣象局于2011年開(kāi)始的“工合(GungHo)”(Globally Uniform Next Generation Highly Optimized)計(jì)劃，以及美國(guó)國(guó)家海洋和大氣管理局(National Oceanic and Atmospheric Administration，NOAA)于2016年開(kāi)始的下一代全球預(yù)報(bào)系統(tǒng)(Next Generation Global Prediction System，NGGPS)計(jì)劃等。

一般而言，要實(shí)現(xiàn)未來(lái)大氣狀態(tài)的氣象要素預(yù)報(bào)，大氣數(shù)值模式需要基于球面計(jì)算網(wǎng)格對(duì)大氣動(dòng)力與熱力學(xué)控制方程進(jìn)行空間數(shù)值離散，繼而采用時(shí)間離散格式向前推進(jìn)模式積分。為此，本文將從數(shù)值算法、準(zhǔn)均勻球面網(wǎng)格和時(shí)間積分方法等3個(gè)方面對(duì)國(guó)內(nèi)外可擴(kuò)展的下一代大氣模式進(jìn)行綜述，希望對(duì)相關(guān)領(lǐng)域同行的研究起到幫助和推動(dòng)作用。需要指出的是，全球模式與區(qū)域模式的差異體現(xiàn)在模式是否有邊界處理，而在數(shù)值算法、球面網(wǎng)格和積分方案方面并無(wú)顯著區(qū)別，因而本文所介紹內(nèi)容對(duì)全球模式與區(qū)域模式均適用。

1 大氣模式數(shù)值方法

自20世紀(jì)20年代RICHARDSON[1]的創(chuàng)始性工作以來(lái)，應(yīng)用在大氣模式中的數(shù)值方法呈現(xiàn)多樣性，如傳統(tǒng)格點(diǎn)模式采用的有限差分方法、有限體積方法、半拉格朗日方法[12-13]和非格點(diǎn)模式(ECWMF)采用的譜方法等。2007年，WILLIAMSON[7]在其綜述“全球大氣模式動(dòng)力框架演變”中詳細(xì)回顧了傳統(tǒng)數(shù)值方法在全球大氣模式中的應(yīng)用進(jìn)展。

近年來(lái)，高性能并行計(jì)算機(jī)的迅速發(fā)展對(duì)應(yīng)用于大氣模擬的數(shù)值計(jì)算方法提出了新的挑戰(zhàn)，即：數(shù)值計(jì)算方法必須能夠充分利用超級(jí)計(jì)算機(jī)所提供的計(jì)算潛力[14-15]?；谶@一背景，過(guò)去十多年中多種適用于大氣數(shù)值模擬并能高效利用超級(jí)計(jì)算機(jī)的高精度數(shù)值方法被應(yīng)用于發(fā)展新一代大氣模式，代表性的方法包括：譜元(spectral element)方法[16-20]、間斷伽遼金(discontinuous Galerkin)方法[21-25]、混合有限元(mixed finite element)方法[26-29]以及多矩約束有限體積(multi-moment constrained finite volume)方法[30-39]。作為有潛力的下一代大氣數(shù)值模式核心基礎(chǔ)算法，這些數(shù)值方法具有高精度、守恒、網(wǎng)格適應(yīng)性強(qiáng)和可擴(kuò)展等特點(diǎn)，國(guó)際先進(jìn)業(yè)務(wù)研究中心開(kāi)始或已經(jīng)采用它們來(lái)發(fā)展新一代大氣數(shù)值模式。與傳統(tǒng)的有限差分和有限體積數(shù)值方法不同，這些數(shù)值方法能夠在局地單元開(kāi)展高階多項(xiàng)式數(shù)值重構(gòu)，因而獲得了高精度和良好可擴(kuò)展特性。然而，上述高精度格式其計(jì)算過(guò)程較為復(fù)雜、數(shù)值計(jì)算魯棒性不高、需要與合適的時(shí)間積分方案相配合等，這些要求給大氣模式開(kāi)發(fā)者采用這類(lèi)高精度數(shù)值方法設(shè)計(jì)下一代數(shù)值模式提出了新的挑戰(zhàn)。

1.1 有限元方法

有限元方法最早發(fā)展于20世紀(jì)40年代初[40]，并在50年代后期取得重要進(jìn)展并得到廣泛應(yīng)用。有限元方法的空間離散格式一般通過(guò)加權(quán)余量法(weighted residual method)導(dǎo)出。在實(shí)際應(yīng)用中，通常采用伽遼金(Galerkin)方法，即：選取離散空間上的基函數(shù)(有時(shí)又稱(chēng)為試探函數(shù)(test function))作為權(quán)重函數(shù)，使之與誤差函數(shù)具有正交性，得到空間離散公式。它的優(yōu)勢(shì)在于能夠靈活地處理復(fù)雜幾何剖分。

有限元方法在固體力學(xué)領(lǐng)域得到了廣泛應(yīng)用[41]，在流體計(jì)算中主要應(yīng)用于低馬赫數(shù)流動(dòng)的模擬[42]。多種有限元方法被應(yīng)用于地球流體力學(xué)數(shù)值模擬，如：基于有限元方法的海洋模式[43-48]。有限元方法應(yīng)用到大氣模擬可追溯到20世紀(jì)60年代[49]，這一領(lǐng)域至20世紀(jì)80年代涌現(xiàn)出更多工作。STANIFORTH[50]對(duì)有限元方法在大氣模擬中的應(yīng)用做過(guò)綜述，強(qiáng)調(diào)有限元模式有望更好地探索、理解和預(yù)報(bào)潛在更復(fù)雜的天氣現(xiàn)象。BéLAND et al.[51]采用有限元方法研究了羅斯貝(Rossby)和重力正規(guī)模態(tài)(normal modes)模擬的準(zhǔn)確性，并且與二階有限差分垂直離散方案做了比較。

相比較于計(jì)算域全局展開(kāi)的譜方法，有限元方法僅在局部單元進(jìn)行有限元函數(shù)級(jí)數(shù)展開(kāi)，沒(méi)有全局并行通信，因而具有良好的計(jì)算效率。然而，與有限差分方法相比，有限元方法空間數(shù)值離散后需要求解代數(shù)方程組，因此它的復(fù)雜程度更高，計(jì)算代價(jià)相對(duì)較大[52]。此外，由于單元內(nèi)未知量多于一個(gè)，有限元方法存在潛在計(jì)算模(computational mode)[53]問(wèn)題，該計(jì)算模是非物理解。特別地，當(dāng)表征小尺度大氣波動(dòng)時(shí)，非線(xiàn)性項(xiàng)的數(shù)值離散容易引起錯(cuò)誤的數(shù)值解。為了抑制計(jì)算模，有限元方法通常需要選擇性地采用耗散機(jī)制[54]，這可能是該方法在大氣模式水平離散中沒(méi)有得到廣泛應(yīng)用的重要原因。相反，采用有限元方法作為垂直離散格式在大氣模式中得到良好應(yīng)用，因?yàn)榇髿鈹?shù)值模式的垂直數(shù)值耗散對(duì)有限元方法抑制計(jì)算模發(fā)揮了關(guān)鍵作用。BURRIDGE et al.[55]在歐洲中期天氣預(yù)報(bào)中心(ECMWF)的預(yù)報(bào)模式中使用有限元方法完成了垂直離散；UNTCH and HORTAL[56]在ECMWF模式中引入垂直有限元平流離散方案，相對(duì)于有限差分模式極大提高了線(xiàn)性重力波位相速度模擬精度，減小了洛倫茲(Lorenz)交叉跳點(diǎn)的計(jì)算模和垂直計(jì)算噪音。

英國(guó)氣象局下一代數(shù)值天氣預(yù)報(bào)模式動(dòng)力框架“工合(GungHo)”研究計(jì)劃將混合有限元方法作為數(shù)值方法的首選。STANIFORTH et al.[28]研究表明，混合有限元離散中速度自由度與氣壓自由度數(shù)量比值小于2時(shí)，會(huì)出現(xiàn)明顯的虛假慣性重力波，而兩者比值大于2時(shí)，則發(fā)現(xiàn)類(lèi)似于六邊形網(wǎng)格上的虛假羅斯貝波[57]。COTTER and SHIPTON[26]針對(duì)線(xiàn)性淺水波方程進(jìn)行了研究，表明混合有限元方法在三角形網(wǎng)格、六邊形網(wǎng)格和正方形網(wǎng)格上通過(guò)選擇適合的基函數(shù)對(duì)(element pair)，調(diào)整氣壓和速度自由度數(shù)量比率，進(jìn)而可以避免虛假的地轉(zhuǎn)模[58]，如：重力慣性模和羅斯貝模。此研究發(fā)現(xiàn)，在三角網(wǎng)格中對(duì)速度離散采用一階Brezzi-Douglas-Fortin-Marini空間而氣壓采用線(xiàn)性不連續(xù)空間，在四邊形網(wǎng)格中速度離散選用k階Raviart-Thomas空間而氣壓則用k階不連續(xù)空間，就可以保證速度自由度和氣壓自由度數(shù)量之比為2∶1，從而避免虛假的地轉(zhuǎn)模，并最小化頻散和耗散誤差。同時(shí)，對(duì)非線(xiàn)性系統(tǒng)連續(xù)方程采用這樣的有限元空間并運(yùn)用間斷伽遼金方法可以保證質(zhì)量守恒。

實(shí)際上，混合有限元方法[59]是跳點(diǎn)網(wǎng)格(如C網(wǎng)格)上有限差分方法的拓展，對(duì)于不同變量選用不同有限元空間函數(shù)(基函數(shù))，規(guī)避了有限差分方法在A格點(diǎn)上的虛假氣壓模，同時(shí)不受數(shù)值網(wǎng)格正交性的約束。目前，這類(lèi)應(yīng)用混合有限元方法在球面淺水波模擬中已獲得良好效果[28-29,60-61]，進(jìn)一步推廣至三維完全可壓縮大氣模式及針對(duì)不同變量的基函數(shù)選擇還在進(jìn)一步研究當(dāng)中[59]。當(dāng)前，MELVIN et al.[62]基于混合有限元離散格式在笛卡爾坐標(biāo)下發(fā)展了三維非靜力動(dòng)力框架，它不僅在模式精度上與現(xiàn)有英國(guó)氣象局業(yè)務(wù)模式ENDGame動(dòng)力框架相當(dāng)，而且具有潛在的眾核高性能計(jì)算與不受模式數(shù)值網(wǎng)格約束的優(yōu)勢(shì)。

1.2 譜元方法

MA[63]最早將譜元方法應(yīng)用于地球流體力學(xué)模擬，建立了適用于復(fù)雜幾何域的譜元淺水模式，研究表明由于譜元方法的耗散和頻散誤差小，模式能夠準(zhǔn)確地模擬海洋鋒面現(xiàn)象。TAYLOR et al.[64]將譜元方法應(yīng)用于球面淺水模式，評(píng)估了譜元方法用于氣候模式的可行性，討論了譜元模式在氣候模擬中的潛在優(yōu)缺點(diǎn)。FOURNIER et al.[19]將譜元方法應(yīng)用到靜力原始方程，建立譜元大氣模式(Spectral Element Atmospheric Model，SEAM)。隨著當(dāng)時(shí)大規(guī)模分布式高性能計(jì)算機(jī)的出現(xiàn)，譜元方法的高可擴(kuò)展性得到了印證。GIRALDO and ROSMOND[18]在水平方向采用譜元方法發(fā)展了一套新的數(shù)值天氣預(yù)報(bào)動(dòng)力框架，測(cè)試了三個(gè)正壓三維標(biāo)準(zhǔn)測(cè)試和四個(gè)斜壓三維標(biāo)準(zhǔn)測(cè)試，如：四波羅斯貝波等；并與四個(gè)天氣預(yù)報(bào)模式和氣候模式相比較，模擬結(jié)果表明譜元大氣模式和譜變換模式精度相當(dāng)，該模式在分布式計(jì)算機(jī)上的并行度可實(shí)現(xiàn)線(xiàn)性可擴(kuò)展。THOMAS and LOFT[65]建立了基于立方球網(wǎng)格的動(dòng)力框架，該框架在水平方向使用高階譜元方法，采用垂直混合氣壓坐標(biāo)和半隱式時(shí)間積分方案。

當(dāng)前，DENNIS et al.[15]發(fā)展的CAM-SE模式、美國(guó)海軍研究院具備譜元可調(diào)的非靜力通用大氣模式(Nonhydrostatic Unified Model for the Atmosphere，NUMA)以及公共地球系統(tǒng)模式(Community Earth System Model，CESM)等均是采用譜元方法發(fā)展的新一代大氣模式。韓國(guó)大氣預(yù)報(bào)系統(tǒng)研究院(Korea Institute of Atmospheric Prediction Systems，KIAPS)在韓國(guó)氣象廳下一代大氣模式中采用高精度譜元方法，具體設(shè)計(jì)方案是在水平方向采用譜元方法離散，在垂直方向采用有限差分離散方法[66-67]。譜元方法具有良好的局地計(jì)算密集型特征，在分布式高性能計(jì)算機(jī)上可獲得高可擴(kuò)展性[15-16]，同時(shí)具有譜收斂和幾何靈活(geometrical flexibility)等優(yōu)勢(shì)，已成為下一代預(yù)報(bào)模式的優(yōu)選數(shù)值計(jì)算方法之一。

1.3 間斷伽遼金方法

間斷伽遼金方法，與前述有限元和譜元方法不同，數(shù)值解和基函數(shù)均局限于單元網(wǎng)格內(nèi)，相鄰網(wǎng)格內(nèi)的近似函數(shù)和權(quán)重函數(shù)在單元界面上的值一般不相等。對(duì)于守恒方程，相鄰單元網(wǎng)格界面的數(shù)值通量一般通過(guò)求解(近似)黎曼(Riemann)問(wèn)題得到，因此有限體積方法中發(fā)展的(近似)黎曼問(wèn)題求解器均可用于間斷伽遼金方法。考慮這一特征，間斷伽遼金方法也被認(rèn)為是一種有限元和有限體積相混合的方法。間斷伽遼金方法和譜元方法一樣也采用不等間距分布的節(jié)點(diǎn)，但其數(shù)值解定義點(diǎn)不包括單元邊界，如采用Legendre-Gauss(LG)點(diǎn)可以獲得非常高的計(jì)算精度。由于間斷伽遼金方法中引入了有限體積的概念，它可以保證數(shù)值守恒性。

REED and HILL[68]最早引入間斷伽遼金方法用于求解靜態(tài)中子線(xiàn)性平流問(wèn)題。LESAINT and RAVIART[69]對(duì)間斷伽遼金方法應(yīng)用于線(xiàn)性平流方程的數(shù)值特性進(jìn)行了分析。CHAVENT and SALZANO[70]首次將間斷伽遼金方法應(yīng)用到求解非線(xiàn)性方程，研究發(fā)現(xiàn)若將一階多項(xiàng)式間斷伽遼金方法空間離散和簡(jiǎn)單向前時(shí)間積分相結(jié)合一般計(jì)算不穩(wěn)定，僅在特定情況下可保證計(jì)算穩(wěn)定性。COCKBURN and SHU[71-72]采用TVB(total variation bounded)限制器、龍格-庫(kù)塔(Runge-Kutta，RK)時(shí)間積分和間斷伽遼金方法求解非線(xiàn)性方程獲得了極大成功。隨后涌現(xiàn)了大量采用間斷伽遼金方法來(lái)求解實(shí)際問(wèn)題的研究工作。SCHWANENBERG and K?NGETER[73]采用龍格-庫(kù)塔間斷伽遼金方法研究了一維淺水方程求解及其在水利上的應(yīng)用，充分展現(xiàn)了其質(zhì)量和動(dòng)量守恒、幾何靈活和緊致性(compactness)好等優(yōu)點(diǎn)。COCKBURN and SHU[74]對(duì)采用龍格-庫(kù)塔時(shí)間積分和間斷伽遼金方法發(fā)展的數(shù)值模型做了較詳細(xì)的綜述。

在大氣模式領(lǐng)域，GIRALDO et al.[75]基于正二十面體網(wǎng)格在笛卡爾坐標(biāo)下發(fā)展了間斷伽遼金淺水波模式，其中空間離散采用nodal類(lèi)型展開(kāi)。NAIR et al.[76-77]基于立方球網(wǎng)格采用modal類(lèi)型展開(kāi)的方式發(fā)展了間斷伽遼金平流模式和淺水波模式，研究結(jié)果表明間斷伽遼金方法的精度與標(biāo)準(zhǔn)譜轉(zhuǎn)換模式相當(dāng)。GIRALDO and RESTELLI[22]將間斷伽遼金方法應(yīng)用于求解通量形式的Navier-Stokes方程，研究了非靜力中尺度大氣現(xiàn)象，數(shù)值試驗(yàn)表明：在熱泡和山脈波試驗(yàn)中間斷伽遼金非靜力模式表現(xiàn)更好；考慮計(jì)算效率，在最壞情況下，間斷伽遼金模式的運(yùn)行速度要比非守恒譜元模式慢近一半，在最好情況下，間斷伽遼金模式的運(yùn)行速度與守恒的譜元模式相當(dāng)。盡管目前采用間斷伽遼金方法建立的大氣模式還未用于實(shí)際業(yè)務(wù)預(yù)報(bào)，然而基于間斷伽遼金方法的研究型大氣模式一直在發(fā)展中，比如美國(guó)國(guó)家大氣研究中心(National Center for Atmospheric Research，NCAR)的HOMME(high-order method modeling environment)靜力模式[24]，美國(guó)海軍研究院的間斷伽遼金全球、區(qū)域一體化NUMA[78]，德國(guó)的DUNE(Distributed and Unified Numerics Environment)大氣模式[79]等。

1.4 多矩約束有限體積方法

有限體積方法是基于物理守恒律發(fā)展的一種數(shù)值方法，它能夠嚴(yán)格保證數(shù)值解的守恒屬性。近年來(lái)人們對(duì)具有良好局地守恒屬性，如緊致模板、高強(qiáng)度浮點(diǎn)計(jì)算和極少數(shù)據(jù)通信的數(shù)值方法產(chǎn)生新的興趣。為此，設(shè)計(jì)能夠適應(yīng)眾核超級(jí)計(jì)算機(jī)的一類(lèi)有限體積方法在非靜力大氣模式[35,80-83]中的應(yīng)用重新受到科學(xué)家們的關(guān)注和研究。

作為對(duì)傳統(tǒng)有限體積法的擴(kuò)展，多矩約束有限體積法是一種數(shù)值求解流體力學(xué)方程組的新方法。多矩方法借鑒了CIP(constrained interpolation profile)方法的基本思想[84-87]，在算法構(gòu)造中同時(shí)采用不同類(lèi)型的矩，如：積分平均值、狀態(tài)變量的點(diǎn)值及其導(dǎo)數(shù)值，不同矩可以采用不同形式的控制方程進(jìn)行預(yù)報(bào)，這些預(yù)報(bào)方程必須保持與原守恒律方程相容。它的構(gòu)建過(guò)程基于明確的物理概念，具有良好可擴(kuò)展性、可達(dá)到高精度、并能保證嚴(yán)格的數(shù)值守恒。與間斷伽遼金和譜元方法相比較，它較好地平衡了數(shù)值模式在計(jì)算精度、魯棒性、算法簡(jiǎn)潔性和計(jì)算效率等各方面的要求，具有較高的實(shí)用價(jià)值。不同于傳統(tǒng)有限體積方法在單元網(wǎng)格內(nèi)定義單個(gè)自由度，多矩約束有限體積方法與譜元和間斷伽遼金方法類(lèi)似，在單元網(wǎng)格內(nèi)定義多個(gè)自由度(或者未知量)，且自由度定義節(jié)點(diǎn)可以自由選擇，包括等距節(jié)點(diǎn)和Legendre-Gauss-Lobatto(LGL)節(jié)點(diǎn)等。該方法通過(guò)施加不同類(lèi)型的約束條件，如：物理量的積分平均值、點(diǎn)值及導(dǎo)數(shù)值，導(dǎo)出自由度的演變方程，最后再通過(guò)龍格-庫(kù)塔時(shí)間積分獲得未知量(或自由度)在新時(shí)刻的預(yù)報(bào)值。由于多矩方法在空間拉格朗日(Lagrange)插值重構(gòu)中引入了多矩約束條件，相比有限元方法中的等距拉格朗日插值具備更好的計(jì)算穩(wěn)定性，能有效抑制高階重構(gòu)函數(shù)在網(wǎng)格單元邊界處的數(shù)值振蕩[32]。XIAO et al.[34]發(fā)展了一套基于多矩約束通量重構(gòu)的通用高精度數(shù)值架構(gòu)，它可以解釋成拉格朗日和埃爾米特(Hermite)混合插值方法，可以從更一般的途徑構(gòu)造高精度數(shù)值格式。同時(shí)研究[88]表明，多矩約束方法在相同自由度下計(jì)算穩(wěn)定的最大庫(kù)朗條件數(shù)(Courant-Friedrichs-Lewy(CFL)數(shù))要大于間斷伽遼金方法。

多矩約束方法已被成功用于發(fā)展地球流體力學(xué)數(shù)值模式。CHEN and XIAO[33]采用多矩約束有限體積方法在立方球面網(wǎng)格上建立了全球淺水波模式。該模式中狀態(tài)變量的點(diǎn)值采用微分形式方程組通過(guò)求解導(dǎo)數(shù)黎曼問(wèn)題來(lái)更新，積分平均值采用通量形式方程組用有限體積公式更新保證數(shù)值守恒性。LI et al.[89]采用球面陰陽(yáng)網(wǎng)格發(fā)展了球面淺水波模式，不同的是狀態(tài)變量點(diǎn)值采用半隱式半拉格朗日方法來(lái)更新。盡管對(duì)狀態(tài)變量點(diǎn)值的更新方法有所不同，但這兩個(gè)淺水模式對(duì)球面大氣波動(dòng)，如：羅斯貝波、過(guò)山氣流波等，均獲得良好的數(shù)值模擬結(jié)果，可與國(guó)際其他先進(jìn)模式相媲美。采用多矩約束有限體積框架，可以構(gòu)造基于局地重構(gòu)的任意高階格式。采用三階、四階和五階多矩約束有限體積方法，基于不同球面準(zhǔn)均勻網(wǎng)格，如：陰陽(yáng)網(wǎng)格、立方球網(wǎng)格、正二十面體三角形網(wǎng)格和六邊形網(wǎng)格，已發(fā)展了各類(lèi)全球淺水模式[36,38-39,90-91]。綜合考慮模式的計(jì)算精度、效率等因素，多矩方法相比其他先進(jìn)數(shù)值方法有其獨(dú)特優(yōu)勢(shì)。

LI et al.[35]首次將多矩約束有限體積方法應(yīng)用于完全可壓縮、非靜力大氣模式。一系列標(biāo)準(zhǔn)數(shù)值試驗(yàn)，如：密度流試驗(yàn)、慣性重力波試驗(yàn)、熱泡試驗(yàn)、Schar復(fù)雜地形試驗(yàn)以及線(xiàn)性靜力/非靜力鐘型山脈波試驗(yàn)[37]等，均表明該模式的計(jì)算精度與當(dāng)前國(guó)際先進(jìn)譜元模式和間斷伽遼金模式相當(dāng)。多矩約束有限體積方法在發(fā)展球面淺水模式和非靜力大氣模式中獲得了良好數(shù)值結(jié)果，有望為發(fā)展我國(guó)下一代高精度、高可擴(kuò)展、守恒的天氣模式甚至氣候模式提供可靠的動(dòng)力學(xué)數(shù)值架構(gòu)。

2 準(zhǔn)均勻球面網(wǎng)格

在大氣數(shù)值模式中代表性的球面網(wǎng)格有：經(jīng)緯度網(wǎng)格、陰陽(yáng)網(wǎng)格、立方球網(wǎng)格和正二十面體網(wǎng)格(圖1)。幾乎所有的業(yè)務(wù)中心均采用過(guò)經(jīng)緯度網(wǎng)格，它具備良好的正交性，方便進(jìn)行差分計(jì)算。然而在超級(jí)高性能計(jì)算環(huán)境中，經(jīng)緯度網(wǎng)格的極區(qū)問(wèn)題及其帶來(lái)的數(shù)據(jù)通信瓶頸變成了制約全球模式性能的一個(gè)關(guān)鍵因素。

圖1 大氣模式中代表性球面網(wǎng)格示意圖(a.經(jīng)緯度網(wǎng)格，b.陰陽(yáng)網(wǎng)格，c.立方球網(wǎng)格，d.正二十面體網(wǎng)格)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]提出了準(zhǔn)均勻、無(wú)數(shù)值奇異的陰陽(yáng)網(wǎng)格，它由兩個(gè)低緯度經(jīng)緯度網(wǎng)格相互扣接而成。陰(或陽(yáng))網(wǎng)格是普通經(jīng)緯度網(wǎng)格一部分，繼承了普通經(jīng)緯度網(wǎng)格特性，且低緯度網(wǎng)格單元面積較均勻，并避免了普通經(jīng)緯度網(wǎng)格數(shù)值計(jì)算奇異性問(wèn)題。球面陰陽(yáng)網(wǎng)格屬于組合網(wǎng)格，在陰陽(yáng)邊界處的數(shù)據(jù)交換需要通過(guò)插值運(yùn)算實(shí)現(xiàn)，為保證質(zhì)量守恒需要特別的考慮[93]。立方球網(wǎng)格由SADOURNY[94]提出，該網(wǎng)格通過(guò)球心投影正立方體表面至球面而成，在其六個(gè)投影面上可生成網(wǎng)格單元面積相當(dāng)均勻的結(jié)構(gòu)化網(wǎng)格。根據(jù)投影方式不同，它可分為心射切面投影(gnomonic)立方球網(wǎng)格和保角(conformal)投影立方球網(wǎng)格[95]。前者生成的球面網(wǎng)格不正交，但擁有解析的坐標(biāo)變換關(guān)系；后者生成的球面網(wǎng)格正交但不存在解析坐標(biāo)變換關(guān)系。SADOURNY et al.[96]提出了正二十面體球面網(wǎng)格，它通過(guò)投影正二十面體至球面得到二十個(gè)球面三角，而在每個(gè)三角內(nèi)又細(xì)分成多個(gè)小三角或六邊形的網(wǎng)格單元，屬于非結(jié)構(gòu)網(wǎng)格。此外，ECWMF從2016年起在業(yè)務(wù)模式中采用立方八面體網(wǎng)格，它是簡(jiǎn)化的高斯網(wǎng)格(reduced Gaussian grid)[97]的一種，即八面體簡(jiǎn)化的高斯網(wǎng)格(octahedral reduced Gaussian grid)，經(jīng)向格點(diǎn)數(shù)(latitude points)保持和簡(jiǎn)化的高斯網(wǎng)格一樣，而相對(duì)于簡(jiǎn)化的高斯網(wǎng)格，其緯向格點(diǎn)數(shù)(longitude points)越往極區(qū)格點(diǎn)數(shù)越少。

網(wǎng)格的幾何屬性在數(shù)值模擬中起到非常重要的作用，它和數(shù)值方法結(jié)合在一起共同影響數(shù)值模擬的準(zhǔn)確性。STANIFORTH and THUBURN[98]對(duì)全球天氣、氣候模式采用的水平網(wǎng)格有詳細(xì)綜述，特別關(guān)注了球面水平網(wǎng)格的計(jì)算模(computational modes)和格點(diǎn)印記(grid imprinting)問(wèn)題。他們認(rèn)為網(wǎng)格幾何屬性會(huì)直接影響數(shù)值模式動(dòng)力框架一些基本的理想性能，如質(zhì)量守恒、準(zhǔn)確描述大氣地轉(zhuǎn)平衡和調(diào)整問(wèn)題、計(jì)算模消失或者被很好地控制住的問(wèn)題、虛假快速傳播的Rossby波問(wèn)題、位勢(shì)和氣壓梯度不應(yīng)產(chǎn)生非物理的渦度源、最小化格點(diǎn)印記等。格點(diǎn)印記是數(shù)值模式底層網(wǎng)格結(jié)構(gòu)潛在誤差的一個(gè)典型現(xiàn)象。數(shù)值試驗(yàn)觀(guān)察發(fā)現(xiàn)，采用準(zhǔn)均勻球面網(wǎng)格的模式其底層網(wǎng)格結(jié)構(gòu)可能會(huì)以噪音或者系統(tǒng)誤差的形式在模擬結(jié)果中凸現(xiàn)出來(lái)，如淺水波試驗(yàn)中立方球網(wǎng)格內(nèi)嵌立方體八個(gè)頂點(diǎn)對(duì)應(yīng)的球面點(diǎn)附近的格點(diǎn)印記[33,38,76-77,99]和陰陽(yáng)網(wǎng)格邊界處格點(diǎn)印記[89,100-102]等。

在當(dāng)前以及未來(lái)超級(jí)高性能計(jì)算條件下，傳統(tǒng)數(shù)值方法，如：有限差分、有限體積方法等，配合經(jīng)緯度網(wǎng)格在數(shù)據(jù)通信問(wèn)題上遇到計(jì)算瓶頸，很可能會(huì)制約它的繼續(xù)發(fā)展。毫無(wú)疑問(wèn)，局地高精度、守恒、高可擴(kuò)展性的數(shù)值方法，如：前述的有限元方法、譜元方法、間斷伽遼金方法以及多矩約束有限體積方法，具有很好的網(wǎng)格靈活性，在準(zhǔn)均勻網(wǎng)格上采用此類(lèi)數(shù)值方法發(fā)展模式動(dòng)力框架是當(dāng)前的一個(gè)熱點(diǎn)研究方向[103]。

3 動(dòng)力框架時(shí)間積分方法

為保證良好的計(jì)算效率，許多全球業(yè)務(wù)天氣預(yù)報(bào)模式采用半隱式半拉格朗日時(shí)間積分方案，比如：平流過(guò)程采用半拉格朗日方法來(lái)計(jì)算，在上游點(diǎn)計(jì)算軌跡不交叉情況下時(shí)間步長(zhǎng)不受限制，而對(duì)于聲波、重力波等快波則采用半隱式方法來(lái)計(jì)算[104]。半隱式模式通過(guò)適當(dāng)?shù)木€(xiàn)性化和空間數(shù)值離散最終得到一個(gè)大規(guī)模、稀疏和對(duì)角化的矩陣，所求解問(wèn)題轉(zhuǎn)換成大型代數(shù)方程，最終采用高效求解器通過(guò)迭代獲得變量的數(shù)值解。

實(shí)際大氣模式網(wǎng)格剖分時(shí)，其水平和垂直網(wǎng)格距之比遠(yuǎn)大于1(十倍甚至于百倍)，為保證計(jì)算效率大氣模式常采用水平顯式-垂直隱式(horizontally explicit-vertically implicit， HEVI)的時(shí)間積分方案，即在水平方向的大氣波動(dòng)項(xiàng)采用顯式積分，而隱式處理垂直傳播波動(dòng)項(xiàng)。相對(duì)于全局半隱式求解方法，它克服了垂直方向的時(shí)間步長(zhǎng)穩(wěn)定性限制，同時(shí)避免了并行計(jì)算中的全局通信。HEVI顯式分裂(split-explicit)方法處理大氣控制方程中不同速度的波動(dòng)時(shí)，將方程中的快變項(xiàng)再分成多個(gè)小時(shí)間步進(jìn)行積分，同時(shí)在垂直方向上用Crank-Nicolson等隱式方法求解[105]，如WRF模式[106]為采用該方法的代表性大氣模式。

大氣模式中的顯式-隱式(implicit-explicit，IMEX)多時(shí)間層方案按照波動(dòng)的快慢分別采用隱式和顯式方案進(jìn)行處理。為了獲得二階或更高精度，一般會(huì)用到多于兩個(gè)時(shí)間層的信息。DURRAN and BLOSSEY[107]用IMEX方法求解靜力原始方程和可壓縮Boussinesq方程。然而，通過(guò)穩(wěn)定性分析可以發(fā)現(xiàn)，多時(shí)間層的IMEX方案會(huì)出現(xiàn)計(jì)算模，因此需要通過(guò)過(guò)濾計(jì)算模以得到滿(mǎn)意的數(shù)值解[108]。相比較而言，IMEX龍格-庫(kù)塔(Runge-Kutta，RK)方案可以基于單時(shí)間層內(nèi)多個(gè)小時(shí)間步獲得良好的穩(wěn)定性和計(jì)算精度(二階及其以上)，即使處理剛性(stiff)物理問(wèn)題[109-110]，在數(shù)學(xué)上表現(xiàn)為算子特征值極大。ULLRICH and JABLONOWSKI[111]采用算子分裂(operator-split)Runge-Kutta-Rosenbrock(RKR)方法積分非靜力大氣模式，在垂直方向采用線(xiàn)性隱式Rosenbrock方法離散，而水平方向采用顯式龍格-庫(kù)塔方法積分，可以獲得三階精度。算子分裂RKR方法實(shí)際上是HEVI方法，因而它的積分穩(wěn)定性最后受限于水平方向網(wǎng)格矩和波動(dòng)速度。值得一提的是，ULLRICH and JABLONOWSKI[111]控制方程在垂直方向所有項(xiàng)均統(tǒng)一采用Rosenbrock方法隱式求解，沒(méi)有區(qū)分波動(dòng)傳播特點(diǎn)，其代價(jià)是會(huì)減慢垂直平流。WELLER et al.[112]采用半隱式或HEVI方法，理論分析了多種龍格-庫(kù)塔IMEX方法求解線(xiàn)性雙曲問(wèn)題的精度和穩(wěn)定性，并提出了一種新的HEVI-split方法，指出該方案可提高重力波模擬的精度和穩(wěn)定性。相對(duì)于半隱式方法，WELLER et al.[112]研究認(rèn)為HEVI-split方法在保持計(jì)算精度上更有優(yōu)勢(shì)。BAO et al.[113]采用基于Strang類(lèi)型的算子分裂HEVI方法積分間斷伽遼金非靜力大氣模式，在水平方向采用魯棒性強(qiáng)保型龍格-庫(kù)塔時(shí)間積分，垂直方向的隱式積分釋放了受聲波約束的時(shí)間步長(zhǎng)，獲得二階時(shí)間收斂精度。

面向未來(lái)并行超級(jí)計(jì)算機(jī)，高階完全顯式時(shí)間積分方案在計(jì)算精度方面無(wú)疑是不錯(cuò)選擇，然而其時(shí)間步長(zhǎng)受限于聲波的快速傳播，當(dāng)前仍然需要在計(jì)算精度和效率上進(jìn)行權(quán)衡。綜合考慮各種因素，在全球模式中采用水平方向顯式求解大氣運(yùn)動(dòng)方程而垂直隱式求解快速波動(dòng)的HEVI方法，并配合上述具有局地高精度特性的數(shù)值算法，在未來(lái)超級(jí)計(jì)算機(jī)平臺(tái)上是一套極具吸引力的地球流體數(shù)值模式架構(gòu)。

4 結(jié)語(yǔ)

世界氣象組織2015年發(fā)布的《Seamless Prediction of the Earth System: From minutes to months》(WMO-No.1156)標(biāo)志著下一代業(yè)務(wù)數(shù)值預(yù)報(bào)的發(fā)展方向在世界范圍內(nèi)形成了共識(shí)。在地球系統(tǒng)科學(xué)框架下，發(fā)展天氣氣候一體化耦合數(shù)值預(yù)報(bào)系統(tǒng)以及相應(yīng)的業(yè)務(wù)預(yù)報(bào)服務(wù)體系成為世界各主要數(shù)值預(yù)報(bào)業(yè)務(wù)中心的研究和應(yīng)用重點(diǎn)。同時(shí)，眾核高性能計(jì)算機(jī)和E級(jí)計(jì)算技術(shù)的發(fā)展也促使下一代業(yè)務(wù)數(shù)值預(yù)報(bào)模式在數(shù)值算法、網(wǎng)格結(jié)構(gòu)等方面做出新的發(fā)展。美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心和德國(guó)氣象局近年來(lái)投入業(yè)務(wù)運(yùn)行的立方球網(wǎng)格有限體積模式(NCEP_GFS-fv3)和基于正二十面體網(wǎng)格的ICON模式是這方面的代表。無(wú)疑，在準(zhǔn)均勻網(wǎng)格上采用局地高精度、守恒、高可擴(kuò)展性的數(shù)值方法是一個(gè)重點(diǎn)和熱點(diǎn)發(fā)展方向。本著拋磚引玉的目的，本文對(duì)近年來(lái)大氣模式中的數(shù)值方法進(jìn)行了系統(tǒng)的分析總結(jié)，希望對(duì)國(guó)內(nèi)同行在下一代數(shù)值預(yù)報(bào)研究中起到幫助和推動(dòng)作用。