蘇勇 沈?qū)W順 彭新東李興良伍湘君張爽陳賢

1南京信息工程大學(xué)大氣科學(xué)學(xué)院，南京210044

2中國氣象科學(xué)研究院災(zāi)害天氣國家重點實驗室，北京100081

3中國氣象局?jǐn)?shù)值預(yù)報中心，北京100081

4廣州空軍氣象中心，廣州510071

1 引言

標(biāo)量平流過程的數(shù)值模擬，一直是數(shù)值模式動力過程發(fā)展的重要內(nèi)容。對于天氣、氣候模式來說，通常需要求解水汽等標(biāo)量的平流輸送方程。因為大氣中水汽的分布具有梯度大、不連續(xù)、多相變等特點，這就給水汽平流過程的精確模擬帶來一定的困難（蘇勇和沈?qū)W順，2009）。一個好的標(biāo)量平流方案一般具備以下幾個特點（Machenhauer et al.,2008）：（1）較高的計算精度；（2）嚴(yán)格的質(zhì)量守恒；（3）正定，沒有明顯的數(shù)值振蕩；（4）保持被平流物理量的形狀；（5）保證計算的穩(wěn)定性；（6）較高的計算效率。

隨著大型計算機(jī)的飛速發(fā)展，數(shù)值模式的分辨率日益提高，這就給平流過程的計算精度提出了更高的要求。對于拉格朗日模式而言，一般情況下，如果直接采用高階精度的插值方案進(jìn)行水汽的平流計算，頻散誤差較大，強(qiáng)梯度區(qū)域的數(shù)值振蕩明顯，不能保證水汽場的空間分布特點和正定、守恒的性質(zhì)（沈?qū)W順等，2011）。

近年來，與計算流體力學(xué)相關(guān)的領(lǐng)域都對標(biāo)量平流過程進(jìn)行了深入的研究。早期的研究主要基于歐拉方法，如 TVD（total variation diminishing）（Thuburn, 1996）方案，基于通量修正思路的FCT（flux-corrected transport）（Book et al., 1975）方案，基于Godunov方法的MUSCL（Monotone Upstreamcentered Schemes for Conservation Laws）（Colella,1985）方案，以及基于分段拋物線函數(shù)的 PPM（Piecewise Parabolic Method）（Colella and Woodward,1984）方案等。其中MUSCL和PPM方案在天氣模式的設(shè)計中得到了廣泛的應(yīng)用和發(fā)展。國內(nèi) Yu（1994）提出了一個兩步保形方案（TSPAS，Twostep Shape Preserving Advection Scheme），有效地改進(jìn)了暴雨中尺度模式中水汽的計算。

歐拉方法容易做到被平流標(biāo)量的守恒性，但也有不易克服的缺點：（1）積分時間步長受到 CFL（Courant-Friedrichs-Lewy）條件的限制；（2）相鄰網(wǎng)格點的插值過程會導(dǎo)致系統(tǒng)的非局地發(fā)展。

20世紀(jì)80年代以后，基于半拉格朗日方法設(shè)計的標(biāo)量平流方案發(fā)展迅速，在天氣模式中得到廣泛應(yīng)用（Robert et al., 1985；Ritchie et al., 1995）。與傳統(tǒng)歐拉方案相比，半拉格朗日方案可取較大得時間步長，積分效率較高，計算精度相當(dāng)，在處理不連續(xù)問題的時候沒有顯著的頻散（Staniforth et al.,1991）。但半拉格朗日平流方案也普遍存在問題，上游點處物理量的內(nèi)插造成總量不能嚴(yán)格守恒（陳嘉濱和季仲貞，2004），在標(biāo)量場的強(qiáng)梯度、不連續(xù)區(qū)域計算精度不高等。

Williamson and Rasch（1989）首先提出半拉格朗日平流方案的保形問題。他們采用三次 Hermite插值或者有理函數(shù)插值，結(jié)合斜率的修正來確保單調(diào)性。Bermejo and Staniforth（1992）設(shè)計出一個混合算法使傳統(tǒng)的半拉格朗日方法轉(zhuǎn)化為準(zhǔn)單調(diào)的半拉格朗日（QMSL）方法，如在平滑的區(qū)域采用三次樣條插值或者三次拉格朗日插值，在不太平滑的區(qū)域采用一階迎風(fēng)格式，或者加大一階迎風(fēng)格式的權(quán)重。這種方法可以抑制強(qiáng)梯度和間斷處的振動，同時在解足夠光滑的區(qū)域保持了傳統(tǒng)半拉格朗日方法的計算精度（沈?qū)W順等，2011），但缺點明顯：（1）在標(biāo)量的強(qiáng)梯度、不連續(xù)區(qū)域，低階插值不能保證計算精度；（2）傳統(tǒng)的拉格朗日算法不能保證標(biāo)量在全場的守恒性。

Xiao et al.（2002）和 Xiao and Peng（2004）從通量形式的平流方程出發(fā)，利用一個網(wǎng)格內(nèi)的積分平均值和兩邊的邊界值，設(shè)計出了PRM（Piecewise Rational Method）方案和 CSLR（Conservative SL with Rational function）方案，其中積分平均值通過通量來更新，邊界值通過半拉格朗日點映射來更新。PRM方案是PPM方案的變形，用有理函數(shù)代替了PPM方案中的拋物線函數(shù),利用有理函數(shù)的保凸性，避免了 PPM 中強(qiáng)制單調(diào)而采用的邊界值調(diào)整。CSLR方案則是基于 CIP-CSL（Constrained Interpolation Profile-Conservative Semi-Lagrangian）（Tanaka et al., 2000; Xiao and Yabe, 2001）方法發(fā)展的一種使物理場精確守恒的半拉格朗日方法，它與 PRM 類似，但在邊界值的計算上，它采用擬合出的曲線中上游點處的值來替換，而 PRM 則采用上游點處周圍格點來插值。這兩個方案無論在標(biāo)量場的平滑處還是間斷處，都可以做到正定保形，并保持較高的精度。雖然 CSLR的算法略簡單于PRM，但CSLR中邊界值是全局變量，在大規(guī)模并行計算中涉及相鄰區(qū)域之間的數(shù)據(jù)交換，程序編制較復(fù)雜，通訊量較大。

中國氣象局自主研發(fā)的新一代全球中期預(yù)報系統(tǒng)GRAPES_GFS（陳德輝和沈?qū)W順，2006；莊世宇等，2005）采用半隱式半拉格朗日的動力框架，標(biāo)量平流過程的計算采用QMSL方案。如前所述，該方案存在不光滑區(qū)域計算精度不高，不能嚴(yán)格守恒的問題。本研究將PRM標(biāo)量平流方案引入GRAPES_GFS中，旨在改進(jìn)模式對水汽平流和分布的模擬，提高降水的預(yù)報效果，提升模式的綜合預(yù)報性能。

第二章簡要介紹PRM標(biāo)量平流方案、GRAPES_GFS中水汽方程以及對極區(qū)采用的技術(shù)處理；第三章通過一維、二維以及模式中的一系列理想試驗對QMSL和PRM兩種方案的性能進(jìn)行對比；第四章通過GRAPES_GFS中的批量預(yù)報試驗，檢驗采用PRM方案之后，對模式預(yù)報效果的改善；第五章給出本研究的結(jié)論。

2 方案簡介

2.1 GRAPES_GFS中的水汽平流方程

GRAPES_GFS大氣運動控制方程組（薛紀(jì)善和陳德輝，2008）中的水汽方程為

其中，q為水汽或其他水物質(zhì)，S(q) 為水物質(zhì)的源匯項，由物理過程部分計算。在不考慮源匯項的情況下，將方程（1）改寫為通量形式：

再將方程（2）推導(dǎo)到球面、地形追隨坐標(biāo)下，得到的水物質(zhì)通量方程為

在上式中散度的推導(dǎo)過程中，暫未考慮地形坡度對散度的影響。u、v分別表示緯向、經(jīng)向速度，表示地形高度追隨坐標(biāo)的高度，表示地形高度追隨坐標(biāo)下的垂直速度，λ和?為經(jīng)度和緯度，a為采用薄層近似之后的地球半徑。

利用分維算法，將方程（3）改寫到三個方向分別進(jìn)行計算：

方程組（4）中三個方程的左邊部分通過一維的平流方程分別計算，右邊部分通過平流之后的散度修正過程計算。

分維算法易于實現(xiàn)，在多維平流及計算中得到廣泛的應(yīng)用，但是當(dāng)流場中在x或y方向有輻合輻散的情況下，該算法會引入虛假的源匯。如果不對這種公式分解過程中產(chǎn)生的差異進(jìn)行處理，那么模式在速度場散度較大的區(qū)域擬合出的水汽場和降水場是不合理的，會有明顯的格點堆積問題，且易引發(fā)積分不穩(wěn)定。Clappier（1998）提出了一種加入修正項的分維算法，在很大程度上緩解了這一問題。本研究采用的即為這種加入修正項的分維算法。

2.2 PRM標(biāo)量平流方案（Xiao and Peng, 2004）

考慮通量形式的一維平流方程：

其中，f為被平流標(biāo)量，t為時間，x為空間坐標(biāo)，u為速度。將上式在控制單元格 [xi?1/2,xi+1/2]內(nèi)積分，經(jīng)過推導(dǎo)可以得到

其中，n為積分時刻，Δx=xi+1/2?xi?1/2，gi±1/2分別代表第i單元格左右邊界的通量。表示第i單元格的積分平均值，

其中，F(xiàn)i(x)表示第i單元格內(nèi)所用算法擬合出的函數(shù)。PRM 算法采用如下的有理函數(shù)擬合計算區(qū)間內(nèi)的標(biāo)量分布形式：

在積分平均值和單元界面值的限制條件下，

fi±1/2分別為第i單元格左、右邊界的值。由（8）、（9）可以求得Ri(x)中參數(shù)的值：

現(xiàn)在，只要再確定了邊界值1/2if±，就可以確定函數(shù)()Fx。PRM 方案利用臨近網(wǎng)格點的積分平均值插值得出邊界值（Colella and Woodward，1984）。

這樣，有理函數(shù)Ri(x)構(gòu)造完成，就可由方程（6）、（7）來計算+1。方程（6）中邊界通量g

i±1/2的表達(dá)式如下:

可以得到：

2.3 極地混合（Peng et al., 2005）

GRAPES_GFS采用球面經(jīng)緯度格點坐標(biāo)，變量在網(wǎng)格上的分布采用 Arakawa-C跳點方案，水物質(zhì)放置在半點上，極點只放置經(jīng)向風(fēng)分量。在球面經(jīng)緯度格點坐標(biāo)中，緯向格距在極點附近變得非常小，如不采取特殊處理會導(dǎo)致積分不穩(wěn)定。因為極點附近水物質(zhì)含量一般不高，且極區(qū)附近四季都以繞極地的緯向風(fēng)為主，經(jīng)向風(fēng)分量不大，所以極點附近水物質(zhì)的分布和變化對高緯度地區(qū)影響不大。由此，引入極地混合技術(shù)，對極點附近小范圍內(nèi)的水物質(zhì)進(jìn)行混合?？紤]到混合區(qū)域與未混合區(qū)域的過渡問題，本研究對傳統(tǒng)的混合技術(shù)做了修正，采取了帶過渡的極地混合。

其中

I為緯向格點總數(shù)，j=1，2，3，4分別對應(yīng)從極點向赤道數(shù)的放置水物質(zhì)的第1、2、3、4個緯圈，q為水物質(zhì)的濃度在格點內(nèi)的積分平均值，avejq為對應(yīng)的第j個緯圈上q的平均值。

采用了考慮過渡的極區(qū)混合技術(shù)之后，PRM 平流方案在極區(qū)的適應(yīng)能力明顯增強(qiáng)，積分穩(wěn)定性也得到了保證。但是當(dāng)極區(qū)經(jīng)向風(fēng)較強(qiáng)時，極區(qū)的混合會帶來一定的偏差，具體可以參考本文GRAPES_GFS中示蹤標(biāo)量穿越極地的理想試驗。

3 理想試驗

將 PRM標(biāo)量平流方案加入 GRAPES_GFS之前，需要通過一系列的理想試驗對其計算精度、頻散性、耗散性等進(jìn)行檢驗，同時與模式中現(xiàn)在采用的QMSL方案進(jìn)行對比。該部分設(shè)計了4組理想試驗：1D方波平流、2D圓錐平流、2D凹槽圓柱旋轉(zhuǎn)和2D變形場試驗。

為了定量檢驗數(shù)值試驗的結(jié)果，定義誤差如下：

通過理想試驗對PRM和QMSL平流方案的性能進(jìn)行了檢驗和對比之后，將 PRM 方案加入GRAPES_GFS的動力框架中，進(jìn)一步通過模式中球面上的理想試驗，對該方案在模式中的合理性，以及在理想情況下的模擬效果進(jìn)行檢驗和對比。

3.1 一維方波平流

一維的區(qū)域內(nèi)有 100個格點，格距為1.0，水平均勻的速度場u=1.0，平流標(biāo)量初始場分布如下：

在不同的時間步長，也就是不同的CFL數(shù)下模擬方波平流 50個格距后的分布形式。這里給出CFL數(shù)取0.5時的模擬結(jié)果。

表1 方波平流試驗計算誤差Table 1 Error of the square wave advection test

從圖1和表1可以看出，QMSL和PRM方案都可以較好地保持方波的形狀，沒有相速度的誤差，都可以完全抑制虛假的上沖和下沖。但相對來說，PRM比QMSL更接近真解，保形效果更好，總誤差、耗散、頻散誤差都比較小。其他CFL條件下的模擬結(jié)果與此相類似。

3.2 凹槽圓柱旋轉(zhuǎn)試驗

二維均勻的區(qū)域在x、y方向都有100個格點，格距為 1.0，區(qū)域內(nèi)為無輻合輻散、繞中心點順時針勻速旋轉(zhuǎn)的速度場，按下式給出：

其中i、j為x、y方向格點數(shù)。初始場如圖2所示，為中心在點（50，50），半徑為25，高為1的帶有凹槽的圓柱，凹槽的范圍為4060,i≤≤2562j≤≤。時間步長取0.5，計算區(qū)域內(nèi)CFL的最大值為0.25，積分12560步，即圓柱順時針旋轉(zhuǎn)10圈。圖2和表2分別給出旋轉(zhuǎn)2圈后和10圈后的結(jié)果以及計算偏差。

圖1 方波平流試驗?zāi)M結(jié)果?？v坐標(biāo)為方波高度; 實線為真解，長虛線為QMSL方案，短虛線為PRM方案Fig.1 The square wave advection test.The ordinate is the height of the square wave.The solid line is the true solution, long-dashed line for QMSL scheme and short-dashed line for PRM scheme

表2 凹槽圓柱旋轉(zhuǎn)試驗計算誤差Table 2 Error of the cut-cylinder rotation test

該試驗主要考察平流方案的保形能力，在沒有輻合輻散的速度場中，初始場的分布形式不應(yīng)該隨著平流過程而發(fā)生變化。從圖2中可以看到，在采用 QMSL方案旋轉(zhuǎn) 2圈之后，已經(jīng)不能夠維持住凹槽圓柱的基本形狀，旋轉(zhuǎn) 10圈之后變形更加明顯。相對來說在采用 PRM 方案的情況下，可以較好的保持初始場中凹槽圓柱的形狀，變形較小。另外從表 2計算誤差的對比中可以看到兩個方案都沒有虛假的下沖，QMSL方案在旋轉(zhuǎn) 10圈之后，區(qū)域內(nèi)的極值被削弱到0.953，而PRM方案還可以達(dá)到 0.999。頻散、耗散誤差的對比同樣類似于前面兩組理想試驗，都是 PRM 方案的誤差明顯較小。

3.3 變形場試驗

本節(jié)采用 QMSL和 PRM方案，重復(fù)了Smolarkiewicz（1982）和 Staniforth and Cote（1987）中的變形場試驗。二維均勻的區(qū)域在x、y方向都有 100個格點，格距為 1.0。計算區(qū)域中采用帶有很強(qiáng)輻合輻散的變形速度場，按照下式給出：

其中i、j為x、y方向格點數(shù)。初始場為中心在點（50，50），底面半徑為15，高為3.87的圓錐。時間步長取0.7，計算區(qū)域內(nèi)CFL的最大值為0.7。參考文獻(xiàn)中給出的試驗結(jié)果，這里也給出分別采用兩種方案的情況下，積分19步、57步、3768步的試驗結(jié)果，如圖3中所示。

該試驗主要考察平流方案在復(fù)雜變形速度場中的適應(yīng)能力，經(jīng)過長時間的積分，數(shù)值解應(yīng)該對稱地分布在中間兩個渦旋的內(nèi)部。這里給出的PRM方案模擬結(jié)果與Staniforth and Cote（1987）中的解析解相近，模擬結(jié)果中體現(xiàn)出了很好的對稱性，但QMSL方案模擬結(jié)果的對稱性稍差，與解析解偏差較大。對比兩個方案的模擬結(jié)果，積分3768步之后，數(shù)值解中的最大值都要小于初始場中的最大值3.87，體現(xiàn)出計算方案的耗散性質(zhì)，所以兩個方案都具有較好的積分穩(wěn)定性，相對來說QMSL方案耗散更大、更穩(wěn)定。

3.4 GRAPES_GFS中的理想平流試驗

在完成上面一系列理想試驗之后，將PRM方案加入GRAPES_GFS的動力框架中，檢驗該方案在模式球面經(jīng)緯度坐標(biāo)中的模擬能力。

圖2 凹槽圓柱旋轉(zhuǎn)試驗?zāi)M結(jié)果。（a）為初始時刻的凹槽圓柱；（b）、（d）分別為QMSL方案旋轉(zhuǎn)2圈和10圈后的結(jié)果；（c）、（e）分別為PRM方案旋轉(zhuǎn)2圈和10圈后的結(jié)果?？v坐標(biāo)為圓柱高度Fig.2 The cut-cylinder rotation test.(a) The initial moment of cut-cylinder; (b), (d) the results of QMSL scheme after 2 laps and 10 laps, respectively; (c),(e) same as (b), (d), but for PRM scheme.The ordinate is the height of the cylinder

圖3 變形場試驗?zāi)M結(jié)果。（a）、（b）、（c）依次為QMSL方案積分19、57、3768步后的模擬結(jié)果；（d）、（e）、（f）依次為PRM方案的對應(yīng)結(jié)果。縱坐標(biāo)為變形后圓錐的高度Fig.3 The deformational flow test.(a), (b), (c) The result of QMSL scheme after 19 steps, 57 steps, and 3768 steps, respectively; (d), (e), (f) same as (a), (b),(c), but for PRM scheme.The ordinate is the height of the cone after deformation

本試驗中，等經(jīng)緯度網(wǎng)格的格距 Δx=Δy= 0.5°,層頂高度ztop= 1.2 km，均勻分為60層。在球面上放置如下定義的示蹤物初始場，如圖4中所示。

其中，λ和?分別為經(jīng)度和緯度，h0= 1.0g/kg，r是球面上距離圓心的大圓距離，水平半徑R取1/3的地球半徑，圓心的垂直高度取為4.5 km，垂直半徑Z取1 km。

驅(qū)動示蹤物進(jìn)行平流的速度場如下設(shè)置：

其中，u0= 38.61 m/s，α為平流方向與赤道的夾角，p0= 1000hPa，運動的尺度H≈ 8.78km。垂直運動的周期為4天，垂直速度的最大值約為4 m/s。初值以及試驗條件的設(shè)置參考 Hundsdorfer and Spee（1995）、Li et al.（,2006）、Jablonowsik et al.（,2008）。

本部分進(jìn)行三組試驗，參考圖4中平流方向的示意，分別對應(yīng)（1）α=0°，示蹤物向西，沿線路1繞地球一周回到初始位置；（2）α=45°，示蹤物向東南，沿線路2繞地球一周回到初始位置；（3）α= 90°，示蹤物向南，沿線路3繞地球兩極一周回到初始位置。時間步長Δt=1200s，積分864步，12天之后示蹤標(biāo)量回到初始位置。

如圖5中所示，在3組試驗中，兩個平流方案的模擬結(jié)果都較為合理，繞地球一周之后，示蹤物的分布形式與初始場可以保持基本一致，沒有位相偏差，中心的極值都有所削弱。但對比之后可以看出，QMSL方案耗散嚴(yán)重，中心極值由1.0耗散至0.6，PRM方案則可以使極值維持在0.8或者0.9；QMSL方案變形嚴(yán)重，x–z剖面的圓形已經(jīng)明顯扭曲，PRM 方案則可以很好地保持住初始場圓形的形狀。

雖然 PRM 方案在前面的理想試驗中表現(xiàn)較好，但是由于采用了極地混合的方案，當(dāng)經(jīng)向風(fēng)分量為主的時候，也就是在90α=°的試驗中，示蹤物過極地時會有明顯的相位以及量值的偏差，如圖 6中水平剖面圖所示。在實際預(yù)報的過程中，當(dāng)過極地氣流較強(qiáng)時，PRM方案也會帶來類似的偏差，這點值得注意。

在沒有源匯項的情況下，標(biāo)量的平流過程應(yīng)該保證區(qū)域內(nèi)總質(zhì)量的守恒。從圖7的對比結(jié)果中可以看到，模式中原來采用的 QMSL方案并不能保證標(biāo)量總質(zhì)量的守恒，在α=0°的試驗中，積分約12天，繞地球一周之后，總質(zhì)量減少了約 1%；在α= 45°的試驗中，總質(zhì)量增加了約 3.5%；在α= 90°的試驗中，兩次穿越極地的過程中總質(zhì)量都有明顯的抖動，誤差最大時可以達(dá)到0.5%。從圖中的變化趨勢可以預(yù)測，如果延長積分時間，在采用QMSL方案的情況下模式中標(biāo)量的總質(zhì)量還會持續(xù)的增大或減小。另一方面，采用歐拉通量形式處理標(biāo)量平流過程的PRM方案則可以做到嚴(yán)格守恒，上面三組試驗中標(biāo)量在全球的總質(zhì)量都不隨積分發(fā)生變化。

通過上述一系列的理想試驗，我們可以看出，PRM 標(biāo)量平流方案相對于模式中原來使用的QMSL方案，精度更高，頻散、耗散誤差更小，正定、保形能力也更強(qiáng)，特別是處理標(biāo)量場的不連續(xù)、強(qiáng)梯度區(qū)域優(yōu)勢明顯；雖然在球面經(jīng)緯度網(wǎng)格系統(tǒng)中需要采用額外的極地混合技術(shù)才能保證穩(wěn)定積分，對過極地標(biāo)量的模擬有一定的影響，但采用歐拉通量形式求解水汽方程的 PRM 方案可以保證積分過程中總質(zhì)量的嚴(yán)格守恒，這一點對全球中期模式，特別是較長時間的積分來說非常重要。

圖4 球面理想試驗初始場及平流方向示意圖Fig.4 Schematic diagram of the initial field and the advection direction of the spherical ideal test

圖5 球面理想試驗，繞地球一周后示蹤物濃度的模擬結(jié)果（x–z剖面，單位：g/kg）。（a）、（b）、（c）依次為QMSL方案在0α=°、45α=°和90α=°的模擬結(jié)果；（d）、（e）、（f）依次為PRM方案的對應(yīng)結(jié)果?？v坐標(biāo)為垂直方向模式層次Fig.5 The simulated concentration of the tracer after a turn around the earth in the spherical ideal test (x–z profile, units: g/kg).(a), (b), (c) The results of QMSL scheme for 0α=°, 45α=°, and 90α=°, respectively; (d), (e), (f) same as (a), (b), (c), but for PRM scheme.The ordinate is the vertical model levels

圖6 PRM方案穿過極地時示蹤物的示意圖Fig.6 The schematic diagram while the tracer crosses the North Pole in PRM scheme

圖7 球面理想試驗中標(biāo)量守恒性的對比。橫坐標(biāo)為繞地球的圈數(shù)，縱坐標(biāo)為每步后全球總質(zhì)量與初值的比例;實線為QMSL方案，虛線為PRM方案Fig.7 Conservation contrast of the spherical ideal test.The abscissa is the number of turns around the earth, the ordinate is the proportion of the global sum after each step to the initial value.Solid line: QMSL scheme;dashed line: PRM scheme

4 批量預(yù)報試驗

上述理想試驗確認(rèn)了 PRM 方案在精度、守恒性及在保形方面的優(yōu)勢，本節(jié)將進(jìn)一步通過在GRAPES_GFS中的實際預(yù)報試驗，檢驗新方案對水物質(zhì)平流、降水預(yù)報效果的改進(jìn)，以及對模式綜合預(yù)報性能的影響。

為此，共進(jìn)行了連續(xù)一個月的8天預(yù)報試驗，時間范圍為2009年7月1日12時（協(xié)調(diào)世界時，下同）至7月31日12時。水平分辨率0.5°，垂直方向36層，時間步長600秒。積分初始場選擇FNL再分析場，不做變分同化過程。試驗所選取的物理過程參數(shù)化方案包括：WSM6微物理方案、RRTMG長短波輻射方案、CoLM（Common Land Model）陸面過程方案、MRF邊界層參數(shù)化方案、SAS積云對流參數(shù)化方案。

下面通過對兩種不同標(biāo)量平流方案情況下濕度場、溫度場、高度場、風(fēng)場，以及水凝物和降水預(yù)報效果的對比，檢驗新標(biāo)量平流方案對模式預(yù)報性能的改善。

從圖8濕度場預(yù)報偏差的對比中可以看出，采用QMSL方案對水物質(zhì)進(jìn)行輸送的時候，模式預(yù)報的中低層濕度場相對于FNL分析場明顯偏濕，如圖8a，72小時預(yù)報的850 hPa濕度場中，亞歐大陸、非洲大陸以及南、北美洲普遍偏濕1～3 g/kg，整個中國東部都是明顯的偏濕區(qū)；在濕度場偏差的緯向平均剖面圖（圖8c）中，北半球低層平均濕度偏大6～8 g/kg，這種偏濕的趨勢隨著對流的運動可以延伸到500 hPa以上。當(dāng)采用PRM方案之后，濕度場正偏差得到明顯的緩解。從圖8b中可以看到，整個亞歐大陸、非洲大陸以及南、北美洲的誤差都得到明顯的控制，圖 8d中模式中層的濕度場正偏差也明顯減小。但是在采用 PRM 方案之后，中低緯地區(qū)太平洋、大西洋以及印度洋上空濕度場的負(fù)偏差有所增大，模擬的大氣明顯偏干，具體原因還有待進(jìn)一步的分析。

標(biāo)量平流方案對模式的影響也不僅限于水物質(zhì)分布和降水的模擬，通過相變加熱、密度變化、各種不穩(wěn)定機(jī)制等復(fù)雜的熱力、動力過程，又影響著模式中溫度場、高度場、風(fēng)場等的預(yù)報。

圖9 批量試驗72小時溫度場預(yù)報平均偏差。單位：K。（a），（b）分別為QMSL方案和PRM方案700 hPa預(yù)報場與FNL的偏差，圖中陰影為偏差，黑線為預(yù)報場；（c），（d）分別為QMSL方案和PRM方案預(yù)報場與FNL偏差的緯向平均垂直剖面Fig.9 Same as Fig.8, but for temperature forecast (units: K)

從圖9a，b溫度場72小時預(yù)報700 hPa偏差的對比中可以看出，采用QMSL方案的情況下，歐亞大陸上空有1～3 K左右的正偏差，采用PRM方案之后上述偏差得到明顯的緩解。從圖9c，d的垂直剖面圖中也可以看到，采用新的PRM方案之后，模式中低層溫度場的預(yù)報偏差有所減小，特別是在北半球中緯度地區(qū)效果明顯。

預(yù)報場相對于分析場的距平相關(guān)系數(shù)（ACC）與均方根偏差（RMSE）可以很好地反映模式的預(yù)報能力，其中，ACC提高表示模式預(yù)報能力增強(qiáng)，RMSE減小表示模式預(yù)報偏差減小。從圖10溫度場北半球的ACC和RMSE可以看出：新的PRM方案對北半球溫度場的預(yù)報有正貢獻(xiàn)，模式的低層、中層、高層 ACC都有所提高；RMSE在中低層有所減小，高層基本保持不變。

進(jìn)一步通過高度場的距平相關(guān)系數(shù)（ACC）與均方根偏差（RMSE）檢驗兩種不同標(biāo)量平流方案對高度場預(yù)報效果的影響。對比圖11和圖10可以看到，采用 PRM 方案之后，相對于給溫度場帶來的變化，對高度場的影響要更加明顯：模式低層、中層、高層高度場距平相關(guān)系數(shù)都明顯提高，均方根偏差都明顯減小。

采用新的 PRM 標(biāo)量平流方案之后，模式低層經(jīng)向風(fēng)場的預(yù)報偏差也有所減小，高原下游地區(qū)低層的南風(fēng)正偏差有所緩解。圖12a中，850 hPa高原東部南風(fēng)明顯偏強(qiáng)，相對于 FNL分析場風(fēng)速偏大5 m/s左右，這種顯著的偏差會對模式造成明顯的影響，如降水帶大范圍北移等。采用 PRM 方案之后，如圖12b所示，高原東部的南風(fēng)偏差減小了約2 m/s，整個華北、東北地區(qū)的風(fēng)場偏差也有所減小。

通過上面的幾組對比可以看到，采用新的PRM標(biāo)量平流方案之后，模式的綜合預(yù)報能力明顯提高，濕度場、溫度場、高度場以及風(fēng)場的預(yù)報偏差都有所減小，最后來檢驗水汽平流過程的改進(jìn)對云的模擬，以及降水預(yù)報效果的影響。

圖10 批量試驗東亞溫度場各層的距平相關(guān)系數(shù)（ACC）與均方根偏差（RMSE）。左邊為ACC，右邊為RMSE。從上至下依次為850 hPa、500 hPa、250 hPa的結(jié)果；實線為QMSL方案；虛線為PRM方案Fig.10 The anomaly correlation coefficient (ACC) and root mean square error(RMSE) of temperature field in East Asia.Solid line with open circles: QMSL scheme; dotted line with solid circles: PRM scheme

圖11 批量試驗東亞高度場各層的距平相關(guān)系數(shù)（ACC）與均方根偏差（RMSE）。左邊為ACC，右邊為RMSE。從上至下依次為850 hPa、500 hPa、250 hPa的結(jié)果；實線為QMSL方案；虛線為PRM方案Fig.11 Same as Fig.10, but for height field

圖12 批量試驗72小時經(jīng)向風(fēng)場預(yù)報平均偏差（單位：m/s）。（a）為QMSL方案850 hPa預(yù)報場與FNL的偏差；（b）為PRM方案850 hPa預(yù)報場與FNL的偏差Fig.12 Average deviation of the 72-h meridional wind forecast at 850 hPa in the bulk test against FNL reanalysis data: (a) QMSL scheme; (b) PRM scheme

前面比較了模式預(yù)報的濕度場相對于 FNL分析場的偏差，但FNL分析場并不能真正的代表水凝物等的實況信息，相對來說，模擬的云帶與衛(wèi)星云圖的對比、預(yù)報降水與實況降水的對比，更能體現(xiàn)出不同水物質(zhì)平流方案在模式中的預(yù)報能力。下面通過對7月21～22日降水過程的個例分析，進(jìn)一步的對比兩種平流方案的預(yù)報性能。

為了與FY2C衛(wèi)星云圖做對比，將模擬的500 hPa以上總水凝物的濃度（云水、雨水、云冰、雪、霰）隨高度進(jìn)行積分（王明歡等，2011），得到模擬云帶的分布，如圖 13a，b所示。將兩種方案模擬的云帶與圖 13c中 FY2C衛(wèi)星云圖對比可以看出，兩種方案在大體上都可以把握住云帶的位置、形狀與量級，相對來說：（1）PRM方案預(yù)報的云帶中心水凝物含量更高，更有利于對強(qiáng)降水過程的模擬；（2）QMSL方案沒有模擬出河套、華北地區(qū)的云帶，PRM方案的模擬結(jié)果與實況更為接近，但云帶位置略有偏北。

另外，在24小時累積降水與實況的對比中可以看出，兩種方案模擬的雨帶在落區(qū)和走勢上差別不大，在強(qiáng)度上略有差異，具體分析可以看到：（1）在強(qiáng)降水中心，PRM方案對降水量級的預(yù)報更加準(zhǔn)確，如江淮流域、印度西北部、黑龍江地區(qū)等；（2）兩種方案模擬的雨帶落區(qū)差別不大，如都漏報了華北北部地區(qū)的降水，長江中下游地區(qū)的雨帶位置都略有偏北等。

對全球中期預(yù)報系統(tǒng)來說，大范圍、長時間系統(tǒng)的模擬能力非常重要。下面對亞洲區(qū)以及全球月平均的24小時累積降水進(jìn)行對比，從圖14中可以看出兩種方案模擬的降水落區(qū)和走勢上差別不大，在局部的強(qiáng)度上略有差異。對比中國區(qū)的月平均降水可以看出，相對于實況，QMSL方案在四川盆地、河套地區(qū)附近模擬出雜亂的大值中心；PRM方案模擬的雨帶比較規(guī)則，與實況更加接近。從這點上可以看出，QMSL方案在地形復(fù)雜、水汽梯度大的區(qū)域模擬效果較差，PRM方案可以更好的處理這一類準(zhǔn)不連續(xù)、大梯度問題。

7月份，赤道輻合帶（ITCZ, intertropical convergence zone）偏向北半球,從西太平洋暖池向東南伸展的南太平洋輻合帶（SPCZ, southern Pacific convergence zone）也比較明顯。從圖14e，f全球的月平均降水中可以看出，兩種方案模擬的降水差異不大，都與實況較為接近；兩種方案都可以把握住ITCZ和SPCZ區(qū)域降水的位置與走勢，只是在局部強(qiáng)度上略有差異。

通過全球 24小時累積降水月平均的緯向平均也可以較為直觀的對比兩種標(biāo)量平流方案下，GRAPES_GFS對各緯度地區(qū)平均降水量的預(yù)報能力。從圖15中可以看出：（1）整體來說，兩種方案模擬的降水沒有大的差異；（2）ITCZ區(qū)域，PRM方案模擬的降水偏強(qiáng)，QMSL方案與實況更為接近；（3）60°S附近和40°N附近，兩種方案模式的降水都較實況偏弱。

下面通過 Threat Score（TS）和 BIAS score（BAIS）評分檢驗，對批量預(yù)報試驗期間各降水量級的降水量以及落區(qū)偏差進(jìn)行定量地檢驗。評分的區(qū)域為東亞區(qū)域（15°～65°N，70°～145°E）。簡要介紹這兩種評分方法（王明歡等，2011）：TS評分主要用來對降水的量級進(jìn)行定量檢驗；bias評分又稱為偏差評分，主要是對降水落區(qū)的偏差進(jìn)行定量檢驗。這兩種評分都是基于累加量級的方法，即分別對大于某一閾值的預(yù)報情況進(jìn)行檢驗，若預(yù)報和實況均為大于此量級的降水即為預(yù)報正確。這兩種評分都是數(shù)值越接近1.0，代表預(yù)報效果越好。

圖13 （a），（b）分別為QMSL和PRM方案2009年7月20日12時24 h預(yù)報的云帶分布（單位：kg/m2）；（c）風(fēng)云2C衛(wèi)星2009年7月21日12時的衛(wèi)星云圖；（d）GPCP資料2009年7月22日00時的24 h累積降水量（單位：mm）；（e），（f）分別為QMSL和PRM方案2009年7月20日12時預(yù)報的12～36 h時段內(nèi)24 h累積降水量（單位：mm）Fig.13 (a), (b) 24-h cloud band forecast at 1200 UTC 20 July 2009 for QMSL and PRM schemes; (c) nephogram of the FY-2C satellite at 1200 UTC 21 July 2009; (d) 24-h accumulated precipitation from GPCP data at 0000 UTC 22 July 2009; (e), (f) the 12–36-h accumulated precipitation at 1200 UTC 20 July 2009 for QMSL and PRM schemes

從圖16東亞區(qū)域降水的檢驗結(jié)果來看，對降水量級的TS檢驗中，對于24 h累積大于10 mm的降水，PRM方案模擬的降水量級更加準(zhǔn)確，但對于降水量在100 mm以上的大暴雨、特大暴雨，模式的預(yù)報能力較差；對落區(qū)偏差的BAIS評分中可以看到，對于各個量級的檢驗，PRM方案的BAIS值都更加接近1.0，模擬的降水落區(qū)都與實況更加接近。

雖然水物質(zhì)平流的模擬對降水非常重要，但它不是決定降水的唯一因素，比如水物質(zhì)之間的相變過程、積云對流、輻射的加熱作用等都在降水的產(chǎn)生過程中發(fā)揮關(guān)鍵的作用。想要更好的模擬降水過程，不僅需要更好的標(biāo)量平流方案，還需要更合理的物理過程參數(shù)化方案。

在這一部分中通過批量實際預(yù)報對 PRM 標(biāo)量平流方案在 GRAPES_GFS中的預(yù)報效果進(jìn)行了檢驗，對比QMSL方案的模擬效果可以看到：（1）PRM方案明顯地減小了中低層濕度場的正偏差；（2）采用 PRM 方案之后，對溫度場、高度場、風(fēng)場的模擬都有正的貢獻(xiàn)；（3）兩種方案模擬的云帶、降水帶沒有大的差異，相對而言，PRM方案對中雨及以上量級的降水模擬的更加準(zhǔn)確，東亞區(qū)的 TS和BAIS評分都更優(yōu)，特別是在地形復(fù)雜、水汽梯度大的區(qū)域模擬的降水帶更加合理。

圖15 批量試驗36～60 h預(yù)報時段內(nèi)24 h累積降水量平均值的緯向平均。黑色實線為GPCP資料；紅色實線為QMSL方案；綠色虛線為PRM方案Fig.15 The zonal average of the averaged 24-h accumulated precipitation in the 36–60-h period in the bulk test.Solid black line: GPCP data; solid red line:QMSL scheme; dashed green line: PRM scheme

圖16 批量試驗東亞區(qū)域降水的TS和BAIS評分。左邊為TS評分，右邊為BAIS評分。從上至下依次為12～36 h、36～60 h和60～84 h的評分結(jié)果。實線為QMSL方案；虛線為PRM方案Fig.16 The Threat Score (TS) and BIAS score (BAIS) of the precipitation in East Asia in the time periods of 12–36 h, 36–60 h, and 60–84 h in the bulk test

5 結(jié)論

PRM方案是Xiao and Peng（2004）利用單元格內(nèi)的積分平均值和兩個邊界值構(gòu)造的基于分段連續(xù)有理函數(shù)的半拉格朗日標(biāo)量平流方案，本研究將其引入GRAPES_GFS，把模式中水汽方程改寫為通量形式后利用其求解，在得到高精度、正定保形數(shù)值解的同時，可以保證積分過程的嚴(yán)格質(zhì)量守恒。在經(jīng)過一系列單機(jī)上的理想試驗、模式動力框架中的理想試驗、批量的實際預(yù)報試驗之后，得到以下幾點主要結(jié)論：

（1）PRM方案相對于QMSL方案，精度更高，頻散、耗散誤差更小，正定、保形能力也更強(qiáng)，特別是處理標(biāo)量場的不連續(xù)、強(qiáng)梯度區(qū)域優(yōu)勢明顯。

（2）GRAPES_GFS中，PRM方案在歐拉格點上求解通量形式的水汽方程，可以保證水物質(zhì)在積分過程中的嚴(yán)格守恒，這一點對全球中期模式，特別是積分時間較長的情況來說非常重要；原來的準(zhǔn)單調(diào)半拉格朗日（QMSL）方案不能保證標(biāo)量的守恒性。

（3）在實際預(yù)報中，相對于QMSL方案，PRM方案明顯地減小了中低層濕度場的正偏差；兩種方案模擬的云帶、降水帶差異不大，相對而言，PRM方案對中雨及以上量級的降水模擬的更加準(zhǔn)確，東亞區(qū)的TS和BAIS評分都更優(yōu)，特別是在地形復(fù)雜、水汽梯度大的區(qū)域模擬的降水帶更加合理。

（4）標(biāo)量平流方案對模式的作用不僅局限于水物質(zhì)和降水的模擬，通過模式中的熱力、動力過程又影響著溫度場、高度場、風(fēng)場等。GRAPES_GFS中采用PRM方案之后，對模式溫度場、高度場、風(fēng)場的模擬都有顯著的正貢獻(xiàn)。

本論文的研究工作取得了一些有意義的結(jié)果，但數(shù)值模式需要靠有限的方程組去求解大氣中復(fù)雜的動力、熱力過程，一個算法在模式中的應(yīng)用也遠(yuǎn)比做理想試驗要復(fù)雜。本研究將 PRM 方案引入GRAPES_GFS之后，也還有一些需要進(jìn)一步解釋和解決的問題：

（1）當(dāng)采用 PRM 方案之后，亞歐大陸、非洲大陸以及南、北美洲地區(qū)低層濕度場正偏差都得到明顯的控制，但中低緯地區(qū)太平洋、大西洋以及印度洋上空濕度場的負(fù)偏差有所增大（圖7b，d），模擬的大氣明顯偏干，具體原因還有待進(jìn)一步的分析。

（2）將PRM方案引入GRAPES_GFS之后，由于經(jīng)緯度網(wǎng)格在極區(qū)的格點輻合問題，需要采用極地混合技術(shù)來保證積分的穩(wěn)定性，這對極地地區(qū)的精確模擬有一定的影響。

（References）

Bermejo R, Staniforth A.1992.The conversion of the semi-Lagrangian advection schemes to quasi-monotone schemes [J].Mon.Wea.Rev., 120:2622–2632.

Book D L, Boris J P, Hain K.1975.Flux-corrected transport II:Generalizations of the method [J].J.Comput.Phys., 18: 248–283.

陳德輝, 沈?qū)W順.2006.新一代數(shù)值預(yù)報系統(tǒng)GRAPES研究進(jìn)展 [J].應(yīng)用氣象學(xué)報, 17: 773–777.Chen D H, Shen X S.2006.Recent progress on GRAPES research and application [J].Journal of Applied Meteorological Science (in Chinese), 17: 773–777.

陳嘉濱, 季仲貞.2004.半隱式半拉格朗日平方守恒格式的構(gòu)造 [J].大氣科學(xué), 28: 527–535.Chen J B, Ji Z Z.2004.A study of complete square-conserving semi-implicit semi-Lagrangian scheme [J].Chinese Journal of Atmospheric Sciences (in Chinese), 28: 527–535.

Clappier A.1998.A correction method for use in multidimensional time-splitting advection algorithms: Application to two- and threedimensional transport [J].Mon.Wea.Rev., 126: 232–242.

Colella P.1985.A direct Eulerian MUSCL scheme for gas dynamics [J].Sci.Stat.Comput., 6: 104–117.

Colella P, Woodward P R.1984.The piecewise parabolic method (PPM) for gas-dynamical simulations [J].J.Comput.Phys., 54: 174–201.

Hundsdorfer W, Spee E J.1995.An efficient horizontal advection scheme for the modeling of global transport of constituents [J].Mon.Wea.Rev.,123: 3554–3564.

Jablonowsik C, Lauritzen P, Nair R, et al.2008.Idealized test cases for the dynamical cores of atmospheric general circulation models: A proposal for the NCAR ASP 2008 summer colloquium.June, 2008, Boulder,Colorado, U.S.A., NCAR.

Li X L, Chen D H, Peng X D, et al.2006.Implementation of the semi-Lagrangian advection scheme on a quasi-uniform overset grid on a sphere [J].Adv.Atmos.Sci., 23: 792–801.

Machenhauer B, Kass E, Lauritzen P H.2008.Finite-Volume Methods in Meteorology [C] // Teman R, Tribbia J, Ciarlet P.Computational Methods for the Atmosphere and the Oceans.New York: Elsevier, 1–120.

Peng X D, Xiao F, Wataru O, et al.2005.Conservative semi-Lagrangian transport on a sphere and the impact on vapor advection in an atmospheric general circulation model [J].Mon.Wea.Rev., 133: 504– 520.

Ritchie H.1995.Implementation of the semi-Lagrangian method in a high resolution version of the ECMWF forecast model [J].Mon.Wea.Rev.,123: 489–514.

Robert A, Yee T L, Ritchie H.1985.A semi-Lagrangian and semi-implicit numerical integration scheme for multilevel atmospheric models [J].Mon.Wea.Rev., 113: 388–394.

沈?qū)W順, 王明歡, 肖峰.2011.GRAPES模式中高精度正定保形物質(zhì)平流方案的研究 I: 理論方案設(shè)計與理想試驗 [J].氣象學(xué)報, 69: 1–15.Shen X S, Wang M H, Xiao F.2011.A study of the high-order accuracy and positive-definite conformal advection scheme in the GRAPES model.Ⅰ: Scientific design and idealized tests [J].Acta Meteorologica Sinica(in Chinese), 69: 1–15.

蘇勇, 沈?qū)W順.2009.相變修正方案在GRAPES模式標(biāo)量平流中的應(yīng)用[J].氣象學(xué)報, 67: 1089–1100.Su Y, Shen X S.2009.Application of amendment to phase-change scheme in scalar advection of GRAPES model [J].Acta Meteorologica Sinica (in Chinese), 67: 1089–1100.

Smolarkiewicz P K.1982.The multi-dimensional Crowley advection scheme [J].Mon.Wea.Rev., 110: 1968–1983.

Staniforth A, C?té J.1991.Semi-Lagrangian integration schemes for atmospheric models—A review [J].Mon.Wea.Rev., 119: 2206–2223.

Staniforth A, Cote J, Pudykiewicz J.1987.Comments on “Smolarkiewicz_s deformational flow” [J].Mon.Wea.Rev., 115: 894–900.

Tanaka R, Nakamura T, Yabe T.2000.Constructing exactly conservative scheme in a non-conservative form [J].Comput.Phys.Commun., 126:232–243.

Thuburn J.1996.TVD schemes, positive schemes, and the universal limiter[J].Mon.Wea.Rev., 125: 1990–1997.

王明歡, 沈?qū)W順, 肖峰.2011.GRAPES模式中高精度正定保形物質(zhì)平流方案的研究II: 連續(xù)實際預(yù)報試驗 [J].氣象學(xué)報, 69: 16–25.Wang M H, Shen X S, Xiao F.2011.A study of the high-order accuracy and positive-definite conformal advection scheme in the GRAPES model.Ⅱ:Continuous actual rainfall prediction experiments [J].Acta Meteorologica Sinica (in Chinese), 69: 16–25.

Williamson D L, Rasch P.1989.Two-dimensional semi-Lagrangian transport with shape preserving interpolation [J].Mon.Wea.Rev., 117:102–129.

Xiao F, Peng X D.2004.A convexity preserving scheme for conservative advection transport [J].J.Comput.Phys., 198: 389–402.

Xiao F, Yabe T.2001.Completely conservative and oscillation-less semi-Lagrangian schemes for advection transportation [J].J.Comput.Phys., 170: 498–522.

Xiao F, Yabe T, Peng X D, et al.2002.Conservative and oscillation-less atmospheric transport schemes based on rational functions [J].J.Geophys.Res., 107: 1–11.

薛紀(jì)善, 陳德輝.2008.數(shù)值預(yù)報系統(tǒng)GRAPES的科學(xué)設(shè)計與應(yīng)用 [M].北京: 科學(xué)出版社, 67–70.Xue J S, Chen D H.2008.Scientific Design and Application on Numerical Prediction System GRAPES (in Chinese)[M].Beijing: Science Press, 67–70.

Yu R C.1994.A two-step shape-preserving advection scheme [J].Adv.Atmos.Sci., 11: 479–490.

莊世宇, 薛紀(jì)善, 朱國富, 等.2005.GRAPES全球三維變分同化系統(tǒng)——基本設(shè)計方案與理想試驗 [J].大氣科學(xué), 29: 872–884.Zhuang S Y, Xue J S, Zhu G F, et al.2005.GRAPES global 3D-Var system—Basic scheme design and single observation test [J].Chinese Journal of Atmospheric Sciences (in Chinese), 29: 872–884.