位溫方程拉格朗日平流計算方案的改進及其在GRAPES_GFS全球模式的應用*

張紅亮 沈學順

國家氣象中心，北京 100081

提 要： GRAPES全球數值預報系統(GRAPES_GFS)個例試驗中，模式高層垂直速度在南極大陸上空會出現較大的計算噪音，進而影響模式的積分穩定性，甚至積分中斷。對GRAPES_GFS的動力框架以及動力、物理耦合方案進行診斷分析，發現位勢溫度在拉格朗日上游點上的垂直插值是垂直速度噪音產生的主要原因。2013年7月實際資料試驗表明，在位溫方程中引入垂直非插值的拉格朗日平流方案，可以減少或消除該計算噪音。使用新的方案后，北半球和熱帶地區模式預報的整體性能得到了明顯提高，8 d預報的全球質量損失較為嚴重的問題也得到了緩解。

引 言

二時間層半隱式半拉格朗日時間積分方案(two-time-level semi-implicit semi-Lagrangian,2TSISL)具有非常好的計算穩定性和計算效率，Robert(1981)提出后就被越來越多地用到了氣象業務模式中，比如:歐洲中期數值預報中心(ECMWF)的數值預報模式IFS(Integrated Forecast System)、英國氣象局的新動力框架模式(New Dynamical Core)以及中國氣象局的GRAPES_GFS模式(Global-Regional Assimilation and PrEdiction System, Global Forecast System) (薛紀善和陳德輝，2008)等。

2009年GRAPES_GFS1.0業務化以來，GRA- PES_GFS在動力框架上做了很多改進，以提高模式的計算穩定性和計算精度。如：地形濾波方案的改進；水物質采用高精度的標量平流方案(piece-wise rational method)；模式高層引入瑞利摩擦；垂直速度方程中引入隱式拖曳項等(沈學順等，2017)；在最新的GRAPES_GFS3.0版本中采用三維參考大氣(蘇勇等，2018)；在時間積分方案上采用預估-修正的2TSISL方案等。

除了上述影響因素，GRAPES_GFS采用的經典2TSISL方案中拉格朗日軌跡的計算精度也是影響模式穩定性和計算精度的一個重要因素。在一個時間步長內，從格點(到達點)出發的空氣質點沿著拉格朗日軌跡反向勻速直線移動所到達的位置就是拉格朗日上游點，經典的2TSISL方案中，質點是以拉格朗日軌跡中間點的速度做勻速直線運動。由于拉格朗日軌跡中間點既不在整時間層上，一般也不在模式格點上，所以質點移動速度需要時間線性外插和空間線性內插得到。Hortal(2002)發現IFS模式采用2TSISL方案后拉格朗日軌跡的計算誤差會造成200 hPa的溫度場和南半球平流層5 hPa急流軸附近的緯向風的計算噪音。他將勻速直線運動假設變成勻加速度直線運動假設，發展了穩定外插的二時間層方案(stable extrapolation two-time-level scheme, SETTLS),一定程度上緩解了這個問題。張旭等(2015)在GRAPES中尺度模式中引入該方案，改進了理想場試驗的計算精度和計算穩定性。由于SETTLS方案本質上仍然是時間外插方案，不能完全消除計算噪音，Wood et al(2014)在英國氣象局的ENDGame中將當前時刻到達點的速度作為質點的移動速度，避免了時間外插帶來的計算噪音，增加了模式的計算穩定性。

GRAPES_GFS預報的垂直速度在南極上空模式高層出現了較大的計算噪音，導致模式積分不穩定，甚至積分中斷的現象。本文從拉格朗日平流的角度對這個噪音的產生進行了深入的診斷分析，得出位溫在拉格朗日上游點上的垂直插值是這個噪音產生的主要原因，進而通過引入位溫垂直非插值拉格朗日平流方案成功地消除了該計算噪音。

1 模式和資料

1.1 2TSISL的GRAPES_GFS簡介

GRAPES_GFS是中國氣象局自主研發的全球區域一體同化預報系統的全球版本。自2009年3月形成準業務化版本GRAPES_GFS1.0以來，針對模式的穩定性和預報精度進行了一系列改造；2014年升級為GRAPES_GFS2.0，模式水平分辨率為0.25°×0.25°，垂直方向不等距分成60層，模式層頂位于36 km；2020年6月將原有的一維參考大氣換成三維參考大氣，經典的2TSISL方案換成預估-修正的2TSISL方案后,GRAPES_GFS3.0正式業務化，模式層頂抬升到60 km，垂直分為87層。

GRAPES_GFS的物理過程包含顯式云微物理方案(Li et al,2018；李喆等，2019)、NSAS對流參數化方案、RRTMG輻射傳輸方案、MRF邊界層過程、CoLM陸面過程方案，以及重力波拖曳等方案(姜曉飛等，2015；萬曉敏等，2017；宮宇等，2018)。

本文試驗是基于GRAPES_GFS2.0進行，其動力框架采用的是經典2TSISL時間積分方案(Temperton et al,2001)。對任意一個預報變量X(緯向風、經向風、垂直速度、無量綱氣壓、水物質和位溫等)，原始預報方程可以用2TSISL展開，得到：

(1)

(2)

線性方程組(1)的求解最后轉為求解關于無量綱擾動氣壓的Helmholtz方程組。 拉格朗日軌跡上游點計算在極區采用Mcdonald and Bates(1989)方法，在中低緯度采用的是Ritchie(1987；1988)方法。拉格朗日軌跡中間點速度為：

(3)

GRAPES_GFS在水平面上采用Arakawa-C網格分布，垂直方向是Charney-Phillips跳點分層，從預報變量緯向風、經向風、無量綱氣壓和位溫(垂直速度等)所在的格點出發需要計算4組拉格朗日軌跡，得到4組拉格朗日上游點位置。拉格朗日上游點上的變量通過32點的三維準三次拉格朗日插值得到。

1.2 資 料

文中采用NCEP/NCAR 提供的水平分辨率為1°×1°，垂直方向26層的FNL(final analysis)資料作為模式的初始場。

2 計算不穩定現象分析

以2011年11月10日12 UTC的FNL資料冷啟個例為例，模式積分4 h后垂直速度在南極上空(>20 000 m)開始出現計算噪音，此時的垂直速度只有0.45 m·s-1，其他預報變量都很光滑。模式積分6 h噪音范圍向北擴大到60°S以南的大部分區域, 位勢溫度和水平風速也相繼出現計算噪音。圖1是GRAPES_GFS積分6 h，第55層(≈22 866 m)模式面上垂直速度分布。垂直速度呈現自西向東正負相間排列分布，最大風速達6 m·s-1，導致模式積分中斷。由于上升和下沉的垂直運動自西向東排列，造成位溫和水平風速自西向東鋸齒狀排列的噪音出現(圖略)。

為甄別計算噪音產生的動力物理原因，關閉全部物理過程后單獨積分動力框架6 h各個變量仍然會出現計算噪音，只是量值減少一半(圖略),因此動力框架的計算誤差是噪音產生的主要原因,物理過程起了放大器作用。

2.1 地形的影響

出現噪音的南極大陸海陸交界處，最大的地形落差有2 000 m(圖2)。GRAPES_GFS采用地形追隨坐標系后，坐標面在此處非常陡峭，因此首先考慮地形對計算噪音的影響。將模式中的地形高度設為零，單獨積分動力框架6 h,第55層模式面上的垂直速度仍然出現了計算噪音(圖3)，最大上升速度為7 m·s-1。所以地形強迫不是計算噪音產生的原因。為了簡化試驗設計，以下試驗中將模式地形設為零。

圖2 研究區域地形高度Fig.2 Topography of research area

圖3 同圖1，但為去掉物理過程和地形Fig.3 Same as Fig.1, but for removing physical packages with topography height being zero

2.2 拉格朗日上游點計算的影響

本節通過敏感性試驗，討論拉格朗日上游點計算對計算噪音的影響。

從拉格朗日到達點出發，沿著i,j,k方向分別計算一個時間步長下質點反向移動的距離得到拉格朗日上游點的位置。GRAPES_GFS 的4套網格點計算4套拉格朗日軌跡并在上游點上進行插值。4套上游點的計算方案和插值方案是一樣的，如果是拉格朗日軌跡計算導致的上游點位置計算誤差帶來的計算噪音，則4套上游點的位置計算都會產生計算噪音。但是上游點上不同變量的插值誤差對模式預報精度的影響是不一樣的。

試驗第一步：討論拉格朗日上游點位置計算誤差對計算噪音的敏感性。分別設置4套拉格朗日平流在i,j,k方向的移動速度為零，在此方向上上游點和到達點處在同一個網格點上，在此方向上也不需要進行空間插值。12組試驗中，只有將垂直速度和位溫所在格點的垂直拉格朗日平流速度設為零時，可以消除計算噪音(圖4)，其他設置下，計算噪音仍然存在。試驗結果表明計算噪音的產生與拉格朗日上游點位置計算無關，垂直速度或位溫在垂直方向的插值誤差是造成計算噪音的主要原因。

試驗第二步：拉格朗日軌跡計算方案不變，垂直速度和位溫在拉格朗日上游點上分別只進行水平插值，不進行垂直插值，以討論垂直速度和位溫的垂直插值誤差對計算噪音的影響。試驗結果表明拉格朗日上游點上的位溫，如果不進行垂直插值而只進行水平插值，計算噪音也會消除(圖5a)；但是，垂直速度不管是否進行垂直插值，計算噪音都會存在。圖5b是垂直速度在上游點上不進行垂直插值的情況，可以看到計算噪音仍然存在，最大上升速度為7 m·s-1。

試驗第三步：在第二步試驗的基礎上，將位溫的拉格朗日上游點在垂直方向上近似到最接近的模式層，拉格朗日上游點上只需要進行水平插值，不做垂直插值。試驗結果表明，垂直速度場是光滑的，計算噪音沒有出現(圖5c)。

以上三步試驗表明，拉格朗日上游點上位溫的垂直插值會造成模式高層垂直速度的計算噪音。

圖4 同圖1，但為位溫所在點的拉格朗日垂直移動速度設為零Fig.4 Same as Fig.1，but for setting the Lagrangian vertical advection velocity on potential temperature points to be zero

圖5 同圖1，但為(a)垂直速度不進行垂直插值,(b)位溫不進行垂直插值, (c)位溫的拉格朗日上游點近似到最近的整數層，位溫不進行垂直插值Fig.5 Same as Fig.1, but for (a) vertical velocity not interpolated vertically; (b) potential temperature not interpolated vertically; (c) approximating the vertical departure point of potential temperature to the nearest model level, with potential temperature not being interpolated vertically

3 垂直非插值拉格朗日平流方案

第2節分析認為，位勢溫度及其相關項在拉格朗日上游點上的垂直插值是模式計算噪音產生的原因。本文據此對位溫方程的拉格朗日垂直平流進行改造,將垂直方向的拉格朗日上游點近似到最接近的模式面,垂直方向不再進行插值。

在沒有源匯項的情況下，位勢溫度θ的熱量方程如下：

(4)

(5)

(6)

rDl=nint(rD)

(7)

式中nint是取整函數。

位溫方程可以重新寫成

(8)

(9)

即：

(10)

4 實際資料試驗

4.1 個例試驗

位溫方程采用垂直非插值拉格朗日平流方案，其他設置和業務模式保持一致，將2011年11月10日12 UTC個例重新積分6 h, 第55層模式面上的垂直速度如圖6所示。由于采用了新的計算方案，各個預報變量上的計算噪音都被消除，最大垂直速度小于0.1 m·s-1。

4.2 批量試驗

以2013年7月1—31日12 UTC 的FNL資料為初值進行8 d預報，模式水平分辨率為0.5°×0.5°，垂直方向60層，積分時間步長為600 s。

31 d試驗中，舊的方案中因為垂直速度出現了計算噪音造成積分中斷的個例共有5次，計算噪音大都發生在南極大陸海陸交界處，或者南北極附近的模式中高層。采用新的垂直平流方案后不穩定現象消失，模式積分8 d正常結束。

圖7為改進預報試驗的預報變量綜合評分卡，分別對東亞、北半球、南半球、赤道四個區域的風場、溫度場、高度場進行檢驗。左列為預報時效1～8 d的距平相關系數，右列為均方根誤差，紅色表示正效果，綠色表示負效果，灰色為中性，箭頭越大效果越顯著?？梢钥吹讲捎眯碌姆桨负螅诒卑肭蚝蜔釒У貐^，模式的綜合預報性能得到明顯提升。在北半球，不論是風場、溫度場或高度場的距平相關系數在1～8 d都有明顯提高，除了250 hPa溫度場的均方根誤差有所變壞，其他預報變量的均方根誤差都有所減少。東亞地區各預報量的距平相關系數與控制試驗相當，溫度場和高度場的均方根誤差有所減少。南半球基本持平。

圖6 同圖1，但為位溫方程采用垂直非插值拉格朗日平流方案Fig.6 Same as Fig.1, but for using the vertical non-interpolated Lagrangian advection scheme in potential temperature equation

圖7 采用新方案后2013年7月模式預報變量評分卡 (紅色：變好,綠色：變差,灰色：相當)Fig.7 Model prediction variable score card of July 2017 (red: better, green: worse, gray: equal)

圖8是2013年7月31 d平均的0～8 d地面氣壓變化。模式起報時平均地面氣壓為982.4 hPa，舊的方案中地面氣壓隨著積分時間線性遞減，平均1 d減少0.2 hPa，經過8 d預報地面氣壓一共減少了1.2 hPa。但是新的方案中8 d預報地面氣壓由982.4 hPa增加到982.7 hPa, 增加了0.3 hPa, 相當于舊方案1.5 d的變化量。

圖8 2013年7月31 d 0～8 d 預報平均地面氣壓變化Fig.8 Averaged surface pressure in 0-8 d over 31 forecasts days starting in July 2013

5 結論與討論

GRAPES_GFS2.0的垂直速度在模式高層出現計算噪音會導致積分不穩定，甚至積分中斷。去掉模式物理過程和模式地形后這個噪音仍然存在，因此認為動力框架計算誤差是計算噪音產生的主要原因。

GRAPES_GFS動力框架采用的是2TSISL方案。拉格朗日軌跡計算包括拉格朗日上游點的位置計算和上游點上的插值。診斷分析表明：拉格朗日上游點位置計算以及其他變量在拉格朗日上游點上的水平和垂直插值，都不會產生垂直速度場的計算噪音，只有位勢溫度在拉格朗日上游點處的垂直插值可以造成垂直速度的計算噪音。在位溫方程中引入垂直非插值拉格朗日平流方案后，將位勢溫度的拉格朗日上游點在垂直方向上近似到最接近的模式層，從而避免了位溫在垂直方向的插值帶來的計算誤差。2013年7月的實際資料試驗表明，新的計算方案可以消除模式中高層的計算噪音，并且可以提高模式在北半球和熱帶地區的預報性能，減少預報的均方根誤差。8 d預報質量損失嚴重的現象得到了明顯的緩解。目前該方案已經在GRAPES_GFS2.0和GRAPES_GFS3.0上業務運行。