PRM標量平流方案在GRAPES全球預報系統中的應用

蘇勇 沈學順 彭新東李興良伍湘君張爽陳賢

1南京信息工程大學大氣科學學院，南京210044

2中國氣象科學研究院災害天氣國家重點實驗室，北京100081

3中國氣象局數值預報中心，北京100081

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

1 引言

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

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

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

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

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

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

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

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

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

2 方案簡介

2.1 GRAPES_GFS中的水汽平流方程

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

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

再將方程（2）推導到球面、地形追隨坐標下，得到的水物質通量方程為

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

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

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

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

2.2 PRM標量平流方案（Xiao and Peng, 2004）

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

其中，f為被平流標量，t為時間，x為空間坐標，u為速度。將上式在控制單元格 [xi?1/2,xi+1/2]內積分，經過推導可以得到

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

其中，Fi(x)表示第i單元格內所用算法擬合出的函數。PRM 算法采用如下的有理函數擬合計算區間內的標量分布形式：

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

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

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

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

i±1/2的表達式如下:

可以得到：

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

GRAPES_GFS采用球面經緯度格點坐標，變量在網格上的分布采用 Arakawa-C跳點方案，水物質放置在半點上，極點只放置經向風分量。在球面經緯度格點坐標中，緯向格距在極點附近變得非常小，如不采取特殊處理會導致積分不穩定。因為極點附近水物質含量一般不高，且極區附近四季都以繞極地的緯向風為主，經向風分量不大，所以極點附近水物質的分布和變化對高緯度地區影響不大。由此，引入極地混合技術，對極點附近小范圍內的水物質進行混合。考慮到混合區域與未混合區域的過渡問題，本研究對傳統的混合技術做了修正，采取了帶過渡的極地混合。

其中

I為緯向格點總數，j=1，2，3，4分別對應從極點向赤道數的放置水物質的第1、2、3、4個緯圈，q為水物質的濃度在格點內的積分平均值，avejq為對應的第j個緯圈上q的平均值。

采用了考慮過渡的極區混合技術之后，PRM 平流方案在極區的適應能力明顯增強，積分穩定性也得到了保證。但是當極區經向風較強時，極區的混合會帶來一定的偏差，具體可以參考本文GRAPES_GFS中示蹤標量穿越極地的理想試驗。

3 理想試驗

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

為了定量檢驗數值試驗的結果，定義誤差如下：

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

3.1 一維方波平流

一維的區域內有 100個格點，格距為1.0，水平均勻的速度場u=1.0，平流標量初始場分布如下：

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

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

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

3.2 凹槽圓柱旋轉試驗

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

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

圖1 方波平流試驗模擬結果。縱坐標為方波高度; 實線為真解，長虛線為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 凹槽圓柱旋轉試驗計算誤差Table 2 Error of the cut-cylinder rotation test

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

3.3 變形場試驗

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

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

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

3.4 GRAPES_GFS中的理想平流試驗

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

圖2 凹槽圓柱旋轉試驗模擬結果。（a）為初始時刻的凹槽圓柱；（b）、（d）分別為QMSL方案旋轉2圈和10圈后的結果；（c）、（e）分別為PRM方案旋轉2圈和10圈后的結果。縱坐標為圓柱高度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 變形場試驗模擬結果。（a）、（b）、（c）依次為QMSL方案積分19、57、3768步后的模擬結果；（d）、（e）、（f）依次為PRM方案的對應結果。縱坐標為變形后圓錐的高度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

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

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

驅動示蹤物進行平流的速度場如下設置：

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

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

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

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

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

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

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

圖5 球面理想試驗，繞地球一周后示蹤物濃度的模擬結果（x–z剖面，單位：g/kg）。（a）、（b）、（c）依次為QMSL方案在0α=°、45α=°和90α=°的模擬結果；（d）、（e）、（f）依次為PRM方案的對應結果。縱坐標為垂直方向模式層次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 球面理想試驗中標量守恒性的對比。橫坐標為繞地球的圈數，縱坐標為每步后全球總質量與初值的比例;實線為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 批量預報試驗

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

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

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

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

標量平流方案對模式的影響也不僅限于水物質分布和降水的模擬，通過相變加熱、密度變化、各種不穩定機制等復雜的熱力、動力過程，又影響著模式中溫度場、高度場、風場等的預報。

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

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

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

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

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

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

圖10 批量試驗東亞溫度場各層的距平相關系數（ACC）與均方根偏差（RMSE）。左邊為ACC，右邊為RMSE。從上至下依次為850 hPa、500 hPa、250 hPa的結果；實線為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 批量試驗東亞高度場各層的距平相關系數（ACC）與均方根偏差（RMSE）。左邊為ACC，右邊為RMSE。從上至下依次為850 hPa、500 hPa、250 hPa的結果；實線為QMSL方案；虛線為PRM方案Fig.11 Same as Fig.10, but for height field

圖12 批量試驗72小時經向風場預報平均偏差（單位：m/s）。（a）為QMSL方案850 hPa預報場與FNL的偏差；（b）為PRM方案850 hPa預報場與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

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

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

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

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

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

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

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

圖13 （a），（b）分別為QMSL和PRM方案2009年7月20日12時24 h預報的云帶分布（單位：kg/m2）；（c）風云2C衛星2009年7月21日12時的衛星云圖；（d）GPCP資料2009年7月22日00時的24 h累積降水量（單位：mm）；（e），（f）分別為QMSL和PRM方案2009年7月20日12時預報的12～36 h時段內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東亞區域降水的檢驗結果來看，對降水量級的TS檢驗中，對于24 h累積大于10 mm的降水，PRM方案模擬的降水量級更加準確，但對于降水量在100 mm以上的大暴雨、特大暴雨，模式的預報能力較差；對落區偏差的BAIS評分中可以看到，對于各個量級的檢驗，PRM方案的BAIS值都更加接近1.0，模擬的降水落區都與實況更加接近。

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

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

圖15 批量試驗36～60 h預報時段內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 批量試驗東亞區域降水的TS和BAIS評分。左邊為TS評分，右邊為BAIS評分。從上至下依次為12～36 h、36～60 h和60～84 h的評分結果。實線為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 結論

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

（1）PRM方案相對于QMSL方案，精度更高，頻散、耗散誤差更小，正定、保形能力也更強，特別是處理標量場的不連續、強梯度區域優勢明顯。

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

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

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

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

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

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

（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.

陳德輝, 沈學順.2006.新一代數值預報系統GRAPES研究進展 [J].應用氣象學報, 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.半隱式半拉格朗日平方守恒格式的構造 [J].大氣科學, 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.

沈學順, 王明歡, 肖峰.2011.GRAPES模式中高精度正定保形物質平流方案的研究 I: 理論方案設計與理想試驗 [J].氣象學報, 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.

蘇勇, 沈學順.2009.相變修正方案在GRAPES模式標量平流中的應用[J].氣象學報, 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.

王明歡, 沈學順, 肖峰.2011.GRAPES模式中高精度正定保形物質平流方案的研究II: 連續實際預報試驗 [J].氣象學報, 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.

薛紀善, 陳德輝.2008.數值預報系統GRAPES的科學設計與應用 [M].北京: 科學出版社, 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.

莊世宇, 薛紀善, 朱國富, 等.2005.GRAPES全球三維變分同化系統——基本設計方案與理想試驗 [J].大氣科學, 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.