我國自主研制的全球/區(qū)域一體化數(shù)值天氣預(yù)報系統(tǒng)GRAPES的應(yīng)用與展望

陳德輝，薛紀(jì)善，沈?qū)W順，孫 健，萬齊林，金之雁，李興良

(1.中國氣象局?jǐn)?shù)值預(yù)報中心，北京 100081;2.中國氣象科學(xué)研究院，北京 100081;3.廣東省氣象局，廣州 510080)

1 前言

氣象災(zāi)害(水災(zāi)、旱災(zāi)、風(fēng)災(zāi)、雪災(zāi)等)與地表災(zāi)害(海嘯、泥石流、潮災(zāi)等)、地質(zhì)構(gòu)造災(zāi)害(地震、火山爆發(fā)等)、生物災(zāi)害(農(nóng)作物和森林的病蟲害、鼠害、傳染病流行等)等其他類型的自然災(zāi)害相比，無論在影響范圍的廣度與深度、影響社會的人數(shù)與活動等方面，還是在造成直接(間接)經(jīng)濟(jì)損失方面，都是最大的，而且一年四季都可能發(fā)生。隨著社會經(jīng)濟(jì)的發(fā)展，氣象災(zāi)害造成的直接經(jīng)濟(jì)損失也在加大。

準(zhǔn)確的天氣預(yù)報是積極防御氣象災(zāi)害、減少人員和財產(chǎn)損失、間接提高生產(chǎn)力的重要手段，而現(xiàn)代天氣預(yù)報則是以數(shù)值天氣預(yù)報技術(shù)作為重要科技支撐。所謂數(shù)值天氣預(yù)報，通俗地說就是計算機(jī)天氣預(yù)報。也就是說，根據(jù)大氣演變的物理規(guī)律(如人們熟知的牛頓運(yùn)動定律、熱力學(xué)定律、質(zhì)量守恒定律等)，應(yīng)用數(shù)值計算的方法，建立數(shù)學(xué)計算模型，編寫成計算機(jī)能自動運(yùn)行計算的軟件，安裝在巨型計算機(jī)上，每天定時輸入氣象觀測數(shù)據(jù)，計算機(jī)就可以在很短的時間內(nèi)制作出未來幾小時、幾天、任意地點(diǎn)的天氣預(yù)報。因此，建立在現(xiàn)代探測技術(shù)、通信技術(shù)、計算機(jī)技術(shù)以及大氣科學(xué)理論基礎(chǔ)上的數(shù)值天氣預(yù)報是解決天氣預(yù)報準(zhǔn)確率持續(xù)提高，滿足社會經(jīng)濟(jì)發(fā)展對防災(zāi)減災(zāi)氣象保障服務(wù)需求的根本科學(xué)途徑。數(shù)值天氣預(yù)報的思想最早是由挪威科學(xué)家Bjecknes于1904提出來的[1]，但是直至第二次世界大戰(zhàn)結(jié)束后，真正意義的“計算機(jī)預(yù)報天氣”才獲得成功[2]。自此以后，數(shù)值天氣預(yù)報準(zhǔn)確率不斷提高，在現(xiàn)代天氣預(yù)報業(yè)務(wù)中發(fā)揮著越來越重要的支撐作用。

以世界上最先進(jìn)的歐洲中期預(yù)報中心(ECMWF)為例，來說明業(yè)務(wù)數(shù)值天氣預(yù)報水平的變化(見表1)。從表1可以看出，從1980年至2010年，500 hPa高度場預(yù)報距平相關(guān)系數(shù)**500 hPa距平相關(guān)系數(shù)是氣象上用于檢驗預(yù)報與實況“相似”程度的一個重要技術(shù)指標(biāo)，數(shù)值越大，兩者的“相似”程度越高，預(yù)報越準(zhǔn)確。一般認(rèn)為距平相關(guān)系數(shù)大于0.6，其預(yù)報結(jié)果是可用的逐年都在提高，比如預(yù)報時效為3 d的距平相關(guān)系數(shù)1980年為84%左右(北半球)，2010年已高達(dá)98%，第5天、第7天、第10天的距平相關(guān)系數(shù)也同樣得到提高，這說明數(shù)值天氣預(yù)報技術(shù)的科學(xué)性及其支撐能力在不斷提升。

表1 ECMWF全球模式的北半球500 hPa高度場預(yù)報距平相關(guān)系數(shù)Table 1 Anomaly correlation coefficient at 500 hPa height field for northern hemisphere at ECMWF

正是基于數(shù)值天氣預(yù)報技術(shù)在現(xiàn)代天氣預(yù)報業(yè)務(wù)中的重要科技支撐作用，世界各國都投入大量的人力物力，布設(shè)地面、雷達(dá)、高空、衛(wèi)星等氣象探測網(wǎng)，購買巨型高性能計算機(jī)，發(fā)展先進(jìn)的數(shù)值天氣預(yù)報系統(tǒng)技術(shù)，為氣象災(zāi)害天氣預(yù)報預(yù)警提供有力支撐。中國氣象局在國家科技研究計劃的支持下，自2001年起，自主研制發(fā)展了我國的新一代全球與區(qū)域一體化數(shù)值天氣預(yù)報系統(tǒng)(GRAPES)[3，4]。

2 GRAPES模式的科學(xué)設(shè)計策略

地球表面的大氣運(yùn)動變化(也即天氣變化)，遵循基本的物理定律:牛頓運(yùn)動定律、熱力學(xué)定律、質(zhì)量守恒定律和大氣狀態(tài)變化關(guān)系[3～5]。

大氣運(yùn)動方程(牛頓運(yùn)動定律):

式(1)～(3)中，u，v，w分別為東西、南北、垂直風(fēng)速分量;r為球半徑(可近似取為地球半徑);f和fφ為科氏力項;λ為經(jīng)度;φ為緯度;t為時間;p為氣壓;δ1和 δ2分別為“0”或“1”參數(shù);Fx如Fu等為次網(wǎng)格物理過程參數(shù)化項。

連續(xù)方程(質(zhì)量守恒定律):

式(4)中，為凈熱源熱匯項;為三維散度;γ=(Cp為定壓比熱，R為干空氣體常數(shù));θ為位溫;Π為氣壓變量。

熱力學(xué)方程(熱力學(xué)定律):

描述溫度T、氣壓p和密度ρ之間變化關(guān)系的大氣狀態(tài)方程為:

式(1)～(8)即為Navies-Stockes方程組。自20世紀(jì)50年代中期以來的 60多年進(jìn)程中[6]，以Navies-Stockes方程組為基礎(chǔ)，各氣象業(yè)務(wù)中心根據(jù)不同的預(yù)報對象需要，研究發(fā)展了多套不同的數(shù)值天氣預(yù)報模式，并同時運(yùn)行中尺度模式、區(qū)域暴雨模式、臺風(fēng)模式、全球預(yù)報模式、短期氣候預(yù)測模式、沙塵暴模式等。這些模式的差異性很大，使得系統(tǒng)運(yùn)行維護(hù)與改進(jìn)開發(fā)的成本高，研究成果的業(yè)務(wù)轉(zhuǎn)化周期長[7，8]。近年來世界各國著手研究建立區(qū)域模式與全球模式統(tǒng)一、天氣與氣候(月、季、年)一體化模式、研究模式與業(yè)務(wù)模式一體化的新一代多尺度統(tǒng)一數(shù)值預(yù)報模式，以取代現(xiàn)有的多套不同尺度、不同預(yù)報對象的模式共存的數(shù)值預(yù)報業(yè)務(wù)體系，大大降低運(yùn)行成本，加速業(yè)務(wù)成果的應(yīng)用轉(zhuǎn)化。目前，加拿大[7，8]、英國[9]、法國[10]等國家的氣象業(yè)務(wù)中心都建立了程度不同的多尺度統(tǒng)一模式系統(tǒng)。

自2001年起，中國氣象局根據(jù)我國氣象數(shù)值天氣預(yù)報業(yè)務(wù)發(fā)展的實際需要，持續(xù)實施了研究發(fā)展我國新一代全球與區(qū)域一體化數(shù)值天氣預(yù)報系統(tǒng)GRAPES的科技攻關(guān)/支撐計劃。

2.1 模式大氣近似假設(shè)

式(1)～(6)方程組是一組不做任何簡化的“完整”的Navies-Stockes方程組，它幾乎包含大氣中各種尺度的運(yùn)動，同時考慮了三維科氏力、球面半徑與球面曲率等影響。GRAPES模式采用淺薄的(r≡a=6 371 km)、“全可壓的”(δ1=1)非靜力平衡大氣的近似假設(shè)(見表2)。

表2 大氣模式方程組的不同簡化假設(shè)方案Table 2 The simple hypothesis methods for equations of the atmospheric model

2.2 模式垂直坐標(biāo)與參考大氣的選擇

2)模式大氣參考廓線。對垂直運(yùn)動方程(3)的量綱分析可知，方程右邊的第一項(氣壓梯度力項)和第二項(重力加速度項)的量級為大項，引入“模式參考大氣:和?表示參考大氣廓線;Π'和θ'表示偏離參考大氣狀態(tài)的擾動量)，消除垂直運(yùn)動方程中滿足靜力平衡的分量，使垂直運(yùn)動方程的差分計算中重力加速度與氣壓梯度力之間由“大量級項”變?yōu)椤皵_動小量級項”，使之降低與方程中其他項的“量級差”，從而有效地提高垂直運(yùn)動方程的計算精度。

2.3 時間離散化方案

GRAPES模式采用半隱式-半拉格朗日時間差分方案，該方案已被世界各大業(yè)務(wù)氣象中心的業(yè)務(wù)數(shù)值天氣預(yù)報模式所廣泛采用，如ECMWF、法國、英國、德國、加拿大、日本、澳大利亞等國家氣象業(yè)務(wù)中心的業(yè)務(wù)數(shù)值預(yù)報模式，其計算公式可簡化為:

式(9)中，ψ 代表 u，v，w，θ，Π 等模式預(yù)報變量;下標(biāo)“D”表示沿拉格朗日軌跡上游出發(fā)點(diǎn)的變量值;右邊表示兩時間層的線性項Lψ和非線性項Nψ權(quán)重平均形式，要求:αε+βε=1。較之歐拉顯式時間差分方案或者時間分裂差分方案，半隱式-半拉格朗日時間差分方案[13]是絕對穩(wěn)定的，理論上可以取很長的時間步長，可大大提高時間差分方案的計算效率，而且具有良好的頻散特性[14]。

對于標(biāo)量(θ和Π)，可以將熱力學(xué)方程和連續(xù)方程直接寫成式(9)的時間離散形式。對于矢量場(u，v和w)的時間離散化，需要采用不同于標(biāo)量場的時間離散化處理。矢量離散法是目前被認(rèn)為較好的方法[15]。采用非線性內(nèi)插三維質(zhì)點(diǎn)位移速度的方法計算拉格朗日質(zhì)點(diǎn)運(yùn)動軌跡，確定上游出發(fā)點(diǎn)。

2.4 空間離散化方案

1)水平網(wǎng)格。作為一個區(qū)域與全球統(tǒng)一的數(shù)值模式，經(jīng)/緯度格點(diǎn)模式是較理想的選擇。一方面，經(jīng)/緯度網(wǎng)格模式的網(wǎng)格設(shè)計簡單，另一方面，可以很方便地設(shè)計有限區(qū)域與全球范圍的嵌套模式(如英國氣象局的統(tǒng)一模式)，也容易與海洋格點(diǎn)模式相耦合。因此，GRAPES模式水平方向采用等經(jīng)-緯度的水平網(wǎng)格(見圖1)和 Arakawa-C格式[16]變量跳點(diǎn)分布設(shè)置(見圖2)，以及二階精度的中央差分格式。

2)垂直離散化。靜力平衡模式由于忽略了垂直運(yùn)動預(yù)報方程(3)，多采用Lorenz的跳層方法設(shè)置模式預(yù)報變量的垂直分布(見圖3)[17]。而對于采用三維運(yùn)動方程的非靜力平衡模式來說，采用非均勻的Charney-Philips跳層方法設(shè)置模式預(yù)報變量的垂直分布[18]，有利于減少從模式動力框架計算到物理過程參數(shù)化計算的溫度場垂直插值計算，減小了垂直氣壓梯度力項的計算誤差，并且可以消除虛假的斜壓不穩(wěn)定[19]。有研究指出[20]，對于處理強(qiáng)渦流與強(qiáng)層流系統(tǒng)、地轉(zhuǎn)適應(yīng)過程等問題，Charney-Philips的跳層[18]設(shè)置優(yōu)于 Lorenz 的跳層[17]設(shè)置。目前，英國氣象局的新一代模式采用Charney-Philips的跳層設(shè)置[18]。因此，GRAPES模式采用Charney-Philips的跳層[18]設(shè)置。

圖1 全球與區(qū)域統(tǒng)一模式的經(jīng)緯度網(wǎng)格系統(tǒng)設(shè)計[3，4]Fig.1 Design of latitude-longitude grid[3，4]for the global-regional unified model

圖2 水平交叉Arakawa-C網(wǎng)格分布[3，4]Fig.2 Arakawa-C staggered grid[3，4]

圖3 變量的Charney-Phillips和Lorenz跳層分布Fig.3 Vertical arrangement of the model prognostic variables with Charney-Phillips’and Lorenz’s methods

2.5 模式程序軟件的設(shè)計策略

為了搭建一個多尺度通用的或者全球/區(qū)域一體化的數(shù)值天氣預(yù)報模式，除了網(wǎng)格系統(tǒng)的選擇、模式大氣的開關(guān)設(shè)置(δ1、δ2)以外，還要從軟件工程構(gòu)造上采取適當(dāng)?shù)牟呗?。參考國外其他國家?shù)值預(yù)報模式程序軟件的設(shè)計策略，尤其是天氣預(yù)報模式(WRF)的程序軟件包基礎(chǔ)架構(gòu)，設(shè)計構(gòu)建一套全新的、統(tǒng)一的GRAPES模式程序軟件包，達(dá)到一套模式程序軟件包，滿足多種研究與業(yè)務(wù)預(yù)報的應(yīng)用需求。

1)標(biāo)準(zhǔn)化。參照其他領(lǐng)域的軟件標(biāo)準(zhǔn)和美國、ECMWF、法國、英國等的數(shù)值天氣預(yù)報系統(tǒng)軟件編程標(biāo)準(zhǔn)，建立了一套既要與國際接軌又要滿足開發(fā)合作與持續(xù)創(chuàng)新的實際需要的、基于FORTRAN-90的氣象數(shù)值預(yù)報模式程序編寫標(biāo)準(zhǔn)。另外，建立一套輔助程序轉(zhuǎn)換實用工具，將已有的FORTRAN-77程序半自動轉(zhuǎn)換為FORTRAN-90程序，以彌補(bǔ)采用手工操作的程序移植方式引起的易出錯、修改不全面等的不足。

2)模塊化。所謂模式程序設(shè)計的模塊化(在FORTRAN-90計算機(jī)語言中，稱為“Module”)就是將數(shù)值天氣預(yù)報系統(tǒng)中的科學(xué)計算內(nèi)容，按照不同的數(shù)值計算方法、計算內(nèi)容、計算順序、不同的物理過程參數(shù)化方案等分為獨(dú)立計算的模塊進(jìn)行編程，達(dá)到根據(jù)不同的需要“插拔式”地選取不同的模塊，組裝成針對不同的預(yù)報對象、不同研究目的的數(shù)值預(yù)報模式，達(dá)到一套模式程序軟件包卻有多種研究與業(yè)務(wù)預(yù)報用途的目的。例如，全球中期數(shù)值天氣預(yù)報模式、有限區(qū)域短期數(shù)值天氣預(yù)報模式、中尺度數(shù)值天氣預(yù)報模式、熱帶臺風(fēng)數(shù)值預(yù)報模式均是由同一套“模塊化”的程序軟件包組裝而成。

3)并行化。多節(jié)點(diǎn)多處理器體系結(jié)構(gòu)的并行計算機(jī)是未來高性能計算機(jī)的主流發(fā)展方向。因此，要使數(shù)值天氣預(yù)報系統(tǒng)在高性能并行計算機(jī)環(huán)境下高效運(yùn)行，對其模式程序作并行化計算處理是不可缺少的重要工作。現(xiàn)有并行化數(shù)值天氣預(yù)報模式程序存在兩方面的問題，一是大部分氣象科學(xué)家編寫的數(shù)值模式，其并行化計算效率不高;二是經(jīng)并行化計算編程專家優(yōu)化的數(shù)值模式，其并行化計算效率高，但模式程序的可讀性較差。而且，由于高性能并行計算機(jī)的快速更換，數(shù)值模式中并行化計算的程序需要作相應(yīng)的大量修改、優(yōu)化。為了更好地發(fā)揮氣象科學(xué)家和并行化計算編程專家的優(yōu)勢作用，并同時考慮到高性能并行計算機(jī)快速更換所帶來的影響，筆者研究團(tuán)隊設(shè)計了將模式程序分為“驅(qū)動層”、“中介層”和“模式層”的“三層結(jié)構(gòu)”的模式程序軟件體系(見圖4)，使并行優(yōu)化與串行模式編程既可獨(dú)立進(jìn)行，又可有機(jī)組合，還可使程序易讀、易寫、易改、易維護(hù)，可移植性強(qiáng)，通用性和靈活性高。如圖4所示，頂層為驅(qū)動層，該層與并行計算機(jī)相關(guān)，底層為模式層，該層程序基本上是串行程序，是大部分氣象科學(xué)家的工作層。在頂層和底層之間為中介層，將驅(qū)動層與模式層有機(jī)地連接起來，在模式程序編譯和運(yùn)行時實現(xiàn)全套模式在負(fù)載較平行的條件下并行化高效計算。

圖4 GRAPES模式程序軟件分層結(jié)構(gòu)Fig.4 The layer-structure of the GRAPES model code software

3 GRAPES模式的應(yīng)用現(xiàn)狀

自2001年以來，在國家“十五”科技攻關(guān)計劃、“十一五”科技支撐計劃的支持下，中國氣象局持續(xù)投入，堅持自主研究發(fā)展了新一代數(shù)值天氣預(yù)報業(yè)務(wù)系統(tǒng)GRAPES。GRAPES模式有兩個系列的應(yīng)用版本:GRAPES_Meso(區(qū)域中尺度模式)和GRAPES_GFS(全球中期預(yù)報模式)。各區(qū)域氣象業(yè)務(wù)中心又在GRAPES_Meso模式的基礎(chǔ)上，發(fā)展針對本區(qū)域天氣特點(diǎn)的業(yè)務(wù)與研究模式，如廣州熱帶海洋氣象研究所(簡稱廣州熱帶所)的GRAPES_TMM(熱帶氣象模式)，上海臺風(fēng)研究所(簡稱上海臺風(fēng)所)的GRAPES_TCM(熱帶氣旋模式)，成都高原氣象研究所(簡稱成都高原所)的GRAPES_TPM(高原氣象模式)等。圖5列舉了GRAPES_Meso和GRAPES_GFS的主要發(fā)展歷程。2001年1月，GRAPES計劃正式啟動。2004年 3月，GRAPES_Meso V1.0(60 km)第一版建立，并在國家氣象中心、廣州熱帶所實時試運(yùn)行應(yīng)用。2005年5月，GRAPES_Meso升級為V1.5(30 km)，在國家氣象中心、廣州熱帶所、上海臺風(fēng)所實時運(yùn)行應(yīng)用。2006年7月，GRAPES_GFS(110 km)全球中期預(yù)報系統(tǒng)建立。2008—2010年，GRAPES_Meso繼續(xù)改進(jìn)升級為 V2.5(150 km)版本，在國家氣象中心、廣州熱帶所、上海臺風(fēng)所、成都高原所業(yè)務(wù)運(yùn)行應(yīng)用。GRAPES_GFS(110 km)全球中期預(yù)報系統(tǒng)正式版本確認(rèn)V1.0(55 km)，并在國家氣象中心業(yè)務(wù)運(yùn)行應(yīng)用。2011年至今，GRAPES區(qū)域模式和全球模式繼續(xù)在改進(jìn)升級、應(yīng)用。表3給出了廣州熱帶所的GRAPES_TMM模式對南海熱帶氣旋(臺風(fēng))路徑預(yù)報誤差的結(jié)果。從表3可以看出，路徑預(yù)報誤差逐年在減小，也即預(yù)報水平在提高，對氣象減災(zāi)防災(zāi)保障服務(wù)起到了很好的科技支撐作用。

圖5 GRAPES模式發(fā)展歷程Fig.5 The historic millstones of GRAPES development

表3 南海熱帶氣旋(臺風(fēng))路徑預(yù)報誤差Table 3 The distance error of tropical cyclones over South China Sea by GRAPES model

GRAPES模式除了在氣象部門的業(yè)務(wù)中心得到應(yīng)用以外，在大氣科學(xué)研究與應(yīng)用培訓(xùn)工作、集合預(yù)報、大氣科學(xué)觀測模擬試驗、沙塵暴預(yù)報、水文洪水預(yù)警預(yù)報、熱帶氣旋預(yù)報、快速更新循環(huán)(RUC)等領(lǐng)域也應(yīng)用十分廣泛，未來還將向氣候模擬、閃電預(yù)報等方向延伸(見圖6)。

圖6 GRAPES模式的主要應(yīng)用領(lǐng)域Fig.6 The main applications of GRAPES

4 GRAPES模式的發(fā)展

GRAPES模式的進(jìn)一步發(fā)展涉及模式初值化、模式動力框架、模式物理過程、專業(yè)應(yīng)用模式拓展、模式不確定性與集合預(yù)報等。這里只就模式動力框架中的垂直坐標(biāo)和網(wǎng)格系統(tǒng)問題進(jìn)行討論。

4.1 模式垂直坐標(biāo)

為了易于處理模式計算域的下邊界條件，氣象數(shù)值預(yù)報模式多采用下邊界光滑的地形追隨坐標(biāo)，其模式面隨地形高低起伏，使得水平氣壓梯度力的計算由原來的一項差分計算變成了兩項差分計算之差，其中一項是與地形坡度相聯(lián)系的項，它會導(dǎo)致水平氣壓梯度力計算的誤差，這種誤差影響從模式低層一直到模式頂層。而且，隨著模式水平分辨率的提高，模式地形坡度越“陡峭”，引起的這類計算誤差越大。對此，有幾種處理方法。第一種是修改氣壓梯度力項的計算方案，例如:回插等壓面法、遞推算法[21]、扣除法[22，23]、基于靜力方程訂正的氣壓回插法[24]等。第二種是 η-階梯地形坐標(biāo)[25]，即用等高面將模式地形切割成“臺階”狀，水平等高面就是模式面。第三種是構(gòu)造新的地形追隨坐標(biāo)面隨高度衰減的函數(shù)，加快“彎曲”的地形追隨坐標(biāo)面隨模式高度增大而逐漸變水平[26]。

第三種方法的優(yōu)點(diǎn)是通過對衰減系數(shù)的調(diào)節(jié)，使得垂直坐標(biāo)在底層受到地形作用的影響逐漸減小，彎曲的坐標(biāo)面逐漸趨緩，到了高層坐標(biāo)面趨于平直，更接近大氣的真實狀況(見圖7)。這一點(diǎn)可以從氣壓梯度力轉(zhuǎn)換計算公式加以理解。

圖7 不同高度的水平模式面Fig.7 Vertical level of model surface

式(10)左邊項為自然高度坐標(biāo)下的氣壓梯度力項，經(jīng)地形追隨坐標(biāo)轉(zhuǎn)換后變成等式右邊的兩項，其中右邊第二項是和地形坡度(b)以及地形坡度雅克比(Jb)有關(guān)的項，bJb隨高度減小，至高層趨于零，這時等式右邊第一項等于左邊項，意味著高層大氣運(yùn)動趨于水平。GRAPES模式將參考第三種方法對垂直坐標(biāo)加以改進(jìn)。

4.2 模式水平網(wǎng)格系統(tǒng)

經(jīng)緯度網(wǎng)格是世界上大多數(shù)國家氣象業(yè)務(wù)中心全球業(yè)務(wù)數(shù)值預(yù)報模式最常用的水平網(wǎng)格系統(tǒng)，其優(yōu)點(diǎn)在于容易構(gòu)造應(yīng)用，而且其正交的經(jīng)緯線可以和南北-東西風(fēng)方向相一致，更易于描述模式大氣的水平運(yùn)動和平流輸送計算。但是，經(jīng)緯度水平網(wǎng)格系統(tǒng)存在致命的缺點(diǎn)。經(jīng)線在南-北兩極輻合(見圖1)，使得靠近極區(qū)的格點(diǎn)距離要大大小于中低赤道地區(qū)的格點(diǎn)距離(如模式水平分辨率取0.1 °×0.1 °時，在赤道附近的格距約為 11.1 km，而在極點(diǎn)附近的格距僅為0.02 km，兩者之比達(dá)573)，并最終在兩極出現(xiàn)奇異點(diǎn)。這些給模式的差分計算帶來了極大的麻煩，解決的辦法如下。

1)精簡格點(diǎn)法。即將靠近極區(qū)處沿緯圈的格點(diǎn)數(shù)減少，而且越靠近極點(diǎn)的緯圈，沿緯圈的格點(diǎn)數(shù)精簡越多。大多數(shù)譜展開模式采用此方法。

2)球冠投影法。即將兩極的區(qū)域(“球冠”)投影到一個均勻網(wǎng)格面上進(jìn)行差分計算，然后再映射回兩極區(qū)域。大多數(shù)格點(diǎn)差分模式采用此方法。

3)設(shè)計一個分布較均勻的水平網(wǎng)格系統(tǒng)，避開傳統(tǒng)的經(jīng)緯度網(wǎng)格經(jīng)線輻合引起的計算不穩(wěn)定的問題，如三角形網(wǎng)格、六角形網(wǎng)格、多邊形網(wǎng)格等。此外，一種靈感可能來自中國古老的道教陰陽八卦理念的新網(wǎng)格系統(tǒng)已被提了出來，這就是所謂的“陰陽”網(wǎng)格系統(tǒng)[27]，或稱疊交網(wǎng)格系統(tǒng)(見圖8c)。

陰陽網(wǎng)格系統(tǒng)由兩個相似的經(jīng)緯度有限區(qū)域網(wǎng)格構(gòu)造成一個完整的全球網(wǎng)格系統(tǒng)，該網(wǎng)格具有格點(diǎn)分布準(zhǔn)均勻和“對稱”、無奇異點(diǎn)的優(yōu)點(diǎn)。該網(wǎng)格的提出為有效解決經(jīng)緯度網(wǎng)格在近極區(qū)格點(diǎn)輻合及極地奇異點(diǎn)計算問題提供了新的思路和途徑。另外，陰陽網(wǎng)格系統(tǒng)的基礎(chǔ)仍然是經(jīng)緯度網(wǎng)格，僅僅是一個坐標(biāo)旋轉(zhuǎn)變換的問題。

式(11)中，x=cosφ·cosλ，y=cosφ·sinλ，z=sinφ;上標(biāo)“i”表示陰網(wǎng)格;上標(biāo)“a”表示陽網(wǎng)格;球坐標(biāo)的計算域為:

式(12)中，δ為設(shè)定的陰陽網(wǎng)格重疊區(qū)。因此，很容易地將經(jīng)緯度系統(tǒng)轉(zhuǎn)換成陰陽網(wǎng)格系統(tǒng)，自然也很容易繼承原來經(jīng)緯度網(wǎng)格的程序軟件。Peng等于2006年成功地實現(xiàn)了在陰陽網(wǎng)格系統(tǒng)上的歐拉平流方案的理論計算[28]。Li等也于同年成功地實現(xiàn)了在陰陽網(wǎng)格系統(tǒng)上的拉格朗日平流計算[29]。GRAPES模式將采用陰陽網(wǎng)格系統(tǒng)以解決傳統(tǒng)經(jīng)緯度網(wǎng)格的極區(qū)計算問題，并嘗試全球-區(qū)域的雙重(或多重)同步雙向嵌套的一體化運(yùn)行策略。

圖8 陰陽網(wǎng)格設(shè)計Fig.8 Design of a Yin-Yang grid

5 結(jié)語

我國科學(xué)家自主研究發(fā)展的GRAPES模式在科學(xué)設(shè)計上采取了單一的動力框架、多種物理過程參數(shù)化方案可選的科學(xué)策略，在軟件工程上采用模塊化、標(biāo)準(zhǔn)化、并行化的軟件程序編碼技術(shù)，使之成為多尺度統(tǒng)一模式，可以用于中尺度、短期區(qū)域、全球中期天氣預(yù)報，區(qū)域與全球氣候預(yù)測，沙塵暴氣溶膠環(huán)境模擬等。水平尺度可從1 km到100 km，時間尺度可從幾小時到2個星期，甚至到月、季、年。這一模式科學(xué)設(shè)計所帶來的優(yōu)點(diǎn)是顯然易見的。

1)使科研開發(fā)與業(yè)務(wù)預(yù)報模式平臺統(tǒng)一，可加快科研成果向業(yè)務(wù)應(yīng)用的轉(zhuǎn)化。

2)可避免一個業(yè)務(wù)中心同時運(yùn)行多套不同構(gòu)架的模式系統(tǒng)的局面，大大降低了業(yè)務(wù)運(yùn)行及維護(hù)成本，包括數(shù)值預(yù)報模式系統(tǒng)在不同體系結(jié)構(gòu)的計算機(jī)上的移植成本。

3)便于更多領(lǐng)域(如數(shù)值模式動力框架、資料同化、物理過程參數(shù)化、軟件工程等)的研發(fā)人員參與，集約化發(fā)展。

作為一個完整的數(shù)值天氣預(yù)報系統(tǒng)，模式積分的初始值生成方案——資料同化系統(tǒng)也是至關(guān)重要的。GRAPES的資料同化系統(tǒng)采用了基于泛函理論的三維/四維變分同化技術(shù)(3D/4D VAR-DA)[30]，以便更有利于同化應(yīng)用更多雷達(dá)、衛(wèi)星遙感觀測資料，構(gòu)造更協(xié)調(diào)的模式積分初值。變分資料同化技術(shù)是目前世界上先進(jìn)的業(yè)務(wù)應(yīng)用技術(shù)。

天氣氣候無縫隙預(yù)報將是未來氣象數(shù)值預(yù)報的發(fā)展方向。針對天氣氣候無縫隙預(yù)報的問題，Shukla[31]指出沒有科學(xué)根據(jù)可以人為地畫出一條邊界來，區(qū)分中尺度預(yù)報、天氣預(yù)報、季節(jié)預(yù)測、ENSO(厄爾尼諾和南方濤動的合稱)預(yù)測，年代際預(yù)測以及氣候變化。但是，從計算機(jī)資源、模式復(fù)雜程度的實際情況考慮，可以建立滿足不同時間尺度需要的天氣氣候預(yù)報預(yù)測系統(tǒng)，而從研究的角度則可以在一個統(tǒng)一的框架下發(fā)展天氣氣候預(yù)報預(yù)測系統(tǒng)。GRAPES模式是多尺度通用的一體化模式，也為未來我國的天氣氣候無縫隙預(yù)報預(yù)測系統(tǒng)的研究發(fā)展奠定了基礎(chǔ)。

[1]Bjecknes V.Das Problem der Wettervorhersage，betrachtet vom Stanpunkt der Mechanik und der Physik[J].Meteor Zeits，1904，21:1-7.

[2]Charney J G，F(xiàn)j?rtoft R，Von Neuman J.Numerical integration of the barotropic vorticity equation[J].Tellus，1950，2:237-254.

[3]陳德輝，薛紀(jì)善，楊學(xué)勝，等.GRAPES新一代全球/區(qū)域多尺度統(tǒng)一數(shù)值預(yù)報模式總體設(shè)計研究[J].科學(xué)通報，2008，53(20):2396-2407.

[4]薛紀(jì)善，陳德輝，沈?qū)W順，等.數(shù)值預(yù)報系統(tǒng)GRAPES的科學(xué)設(shè)計與應(yīng)用[M].北京:科學(xué)出版社，2008.

[5]Haltiner G J，Williams R T.Numerical Prediction and Dynamic Meteorology[M].2nd ed.New York:John Wiley＆ Sons Press，1980:1-3.

[6]Kalnay E.Atmospheric Modeling，Data Assimilation and Predictability[M].Cambridge:Cambridge University Press，2003:4-12.

[7]Cote J，Gravel S，Methot A，et al.The operational CMC-MRB Global Environmental Multiscale(GEM)model.Part I:Design considerations and formulation[J].Monthly Weather Review，1998，126:1373-1395.

[8]Cote J，Gravel S，Methot A，et al.The operational CMC-MRB global environmental multiscale(GEM)model.Part II:Results[J].Monthly Weather Review，1998，126:1397-1418.

[9]Cullen M J P.The unified forecast/climate model[J].Meteorological Magazine，1993，122:81-94.

[10]Bubnova R，Gello G，Bernard P，et al.Integration of the fully elastic equations cast in hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/ALADIN NWP system[J].Monthly Weather Review，1995，123:515-535.

[11]Staniforth A，Cote J J.Semi-lagrangian integration schemes for atmospheric models review [J].Monthly Weather Review，1991，119:2206-2223.

[12]Gal-Chen T，Sommerville R C J.On the use of a coordinate transformation for the solution of the Navies-Stokes equations[J].Journal of Computing Physics，1975，17:209-228.

[13]Robert A.A semi lagrangian and semi implicit numerical integration scheme for the primitive meteorological equations[J].Journal of the Meteorological Society of Japan，1982，60:319-325.

[14]Hortal M.Aspects of the numeric of the ECMWF model[C]//Proceedings of ECMWWF Workshop on Recent Development in Numerical Methods for Atmospheric Modelling.U.K.:Shinfield Park，1999:127-143.

[15]Qian J H，Semazzi F H M，Scroggs J S.A global nonhydrostatic semi-Lagrangian atmospheric model with orography[J].Monthly Weather Review，1998，126:747-770.

[16]Arakawa A，Lamb V R.Computational design of the basic dynamical processes of the UCLA general circulation model[M]//Methods of Computational Physics，General Circulation Models of the Atmosphere.New York:Academic Press，1977:174-265.

[17]Lorenz E N.The predictability of a flow which possesses many scales of motion[J].Tellus，1969，21:289-307.

[18]Charney J G，Phillips N A.Numerical integration of the quasigeostrophic equations for barotronic and simple baroclinic flows[J].Journal of Meteorology，1953，10:71-99.

[19]Arakawa A，Konor C S.Vertical differencing of the primitive equation based on the Charney-Phillips grid in Hybrid σ-p vertical coordinates[J].Monthly Weather Review，1996，124:511-528.

[20]Cullen M J P.The use of dynamical knowledge of the atmosphere to improve NWWP models[C]//Proceedings of ECMWF Workshop on Recent Development in Numerical Methods for Atmospheric Modelling.U.K.:Shinfield Park，1999:418-441.

[21]楊曉娟，錢永甫.p-σ坐標(biāo)系數(shù)值模式中氣壓梯度力的遞推算法[J].大氣科學(xué)，2003，27(2):171-181.

[22]錢永甫，王云峰.數(shù)值模式中氣壓梯度力的算法試驗[J].氣象學(xué)報，1991，49(3):538-547.

[23]周天軍，錢永甫.地形區(qū)兩種典型氣壓梯度力計算方法的比較[J].應(yīng)用氣象學(xué)報，1996，7(3):377-380.

[24]胡江林，王盤興.地形跟隨坐標(biāo)下的中尺度模式氣壓梯度力計算誤差分析及其改進(jìn)方案[J].大氣科學(xué)，2007，31(1):110-118.

[25]Mesinger F.A blocking technique for representation of mountains in atmospheric models[J].Rivista di Meteorologia Aeronautica，1984，44:195-202.

[26]Schar C，Leuenberger D，F(xiàn)uhrer O，et al.A New terrain-following vertical coordinate formulation for atmospheric prediction models[J].Monthly Weather Review，2002，130:2459-2480.

[27]Kageyama A，Sato T.The“Yin-Yang grid”:an overset grid in spherical geometry[J].Geochemistry，Geophysics，Geosystems，2004，5:1-15.

[28]Peng X， Xiao F，Takahashi K.Conservative constraint for a quasi-uniform overset grid on sphere[J].Quarterly Journal of the Royal Meteorological Society，2006，132:979-996.

[29]Li Xingliang，Chen Dehui，Peng Xindong，et al.Implementation of the semi-Lagrangian advection scheme on a quasi-uniform overset grid on a sphere [J].Advances in Atmospheric Sciences，2006，23(5):792-801.

[30]薛紀(jì)善，莊世宇，朱國富，等.新一代全球/區(qū)域多尺度通用變分同化系統(tǒng)研究[J].科學(xué)通報，2008，53(20):2407-2417.

[31]Shukla J.Seamless prediction of weather and climate:a new paradigm for modeling and prediction research[C]//U.S.National Oceanic and Atmospheric Administration“Climate Test Bed Joint Seminar Series”，NCEP.Maryland:Camp Springs，2009:10.