位溫方程拉格朗日平流計(jì)算方案的改進(jìn)及其在GRAPES_GFS全球模式的應(yīng)用*

張紅亮 沈?qū)W順

國家氣象中心，北京 100081

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

引 言

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

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

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

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

1 模式和資料

1.1 2TSISL的GRAPES_GFS簡(jiǎn)介

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

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

本文試驗(yàn)是基于GRAPES_GFS2.0進(jìn)行，其動(dòng)力框架采用的是經(jīng)典2TSISL時(shí)間積分方案(Temperton et al,2001)。對(duì)任意一個(gè)預(yù)報(bào)變量X(緯向風(fēng)、經(jīng)向風(fēng)、垂直速度、無量綱氣壓、水物質(zhì)和位溫等)，原始預(yù)報(bào)方程可以用2TSISL展開，得到：

(1)

(2)

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

(3)

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

1.2 資 料

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

2 計(jì)算不穩(wěn)定現(xiàn)象分析

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

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

2.1 地形的影響

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

圖2 研究區(qū)域地形高度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 拉格朗日上游點(diǎn)計(jì)算的影響

本節(jié)通過敏感性試驗(yàn)，討論拉格朗日上游點(diǎn)計(jì)算對(duì)計(jì)算噪音的影響。

從拉格朗日到達(dá)點(diǎn)出發(fā)，沿著i,j,k方向分別計(jì)算一個(gè)時(shí)間步長(zhǎng)下質(zhì)點(diǎn)反向移動(dòng)的距離得到拉格朗日上游點(diǎn)的位置。GRAPES_GFS 的4套網(wǎng)格點(diǎn)計(jì)算4套拉格朗日軌跡并在上游點(diǎn)上進(jìn)行插值。4套上游點(diǎn)的計(jì)算方案和插值方案是一樣的，如果是拉格朗日軌跡計(jì)算導(dǎo)致的上游點(diǎn)位置計(jì)算誤差帶來的計(jì)算噪音，則4套上游點(diǎn)的位置計(jì)算都會(huì)產(chǎn)生計(jì)算噪音。但是上游點(diǎn)上不同變量的插值誤差對(duì)模式預(yù)報(bào)精度的影響是不一樣的。

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

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

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

以上三步試驗(yàn)表明，拉格朗日上游點(diǎn)上位溫的垂直插值會(huì)造成模式高層垂直速度的計(jì)算噪音。

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

圖5 同圖1，但為(a)垂直速度不進(jìn)行垂直插值,(b)位溫不進(jìn)行垂直插值, (c)位溫的拉格朗日上游點(diǎn)近似到最近的整數(shù)層，位溫不進(jìn)行垂直插值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節(jié)分析認(rèn)為，位勢(shì)溫度及其相關(guān)項(xiàng)在拉格朗日上游點(diǎn)上的垂直插值是模式計(jì)算噪音產(chǎn)生的原因。本文據(jù)此對(duì)位溫方程的拉格朗日垂直平流進(jìn)行改造,將垂直方向的拉格朗日上游點(diǎn)近似到最接近的模式面,垂直方向不再進(jìn)行插值。

在沒有源匯項(xiàng)的情況下，位勢(shì)溫度θ的熱量方程如下：

(4)

(5)

(6)

rDl=nint(rD)

(7)

式中nint是取整函數(shù)。

位溫方程可以重新寫成

(8)

(9)

即：

(10)

4 實(shí)際資料試驗(yàn)

4.1 個(gè)例試驗(yàn)

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

4.2 批量試驗(yàn)

以2013年7月1—31日12 UTC 的FNL資料為初值進(jìn)行8 d預(yù)報(bào)，模式水平分辨率為0.5°×0.5°，垂直方向60層，積分時(shí)間步長(zhǎng)為600 s。

31 d試驗(yàn)中，舊的方案中因?yàn)榇怪彼俣瘸霈F(xiàn)了計(jì)算噪音造成積分中斷的個(gè)例共有5次，計(jì)算噪音大都發(fā)生在南極大陸海陸交界處，或者南北極附近的模式中高層。采用新的垂直平流方案后不穩(wěn)定現(xiàn)象消失，模式積分8 d正常結(jié)束。

圖7為改進(jìn)預(yù)報(bào)試驗(yàn)的預(yù)報(bào)變量綜合評(píng)分卡，分別對(duì)東亞、北半球、南半球、赤道四個(gè)區(qū)域的風(fēng)場(chǎng)、溫度場(chǎng)、高度場(chǎng)進(jìn)行檢驗(yàn)。左列為預(yù)報(bào)時(shí)效1～8 d的距平相關(guān)系數(shù)，右列為均方根誤差，紅色表示正效果，綠色表示負(fù)效果，灰色為中性，箭頭越大效果越顯著。可以看到采用新的方案后，在北半球和熱帶地區(qū)，模式的綜合預(yù)報(bào)性能得到明顯提升。在北半球，不論是風(fēng)場(chǎng)、溫度場(chǎng)或高度場(chǎng)的距平相關(guān)系數(shù)在1～8 d都有明顯提高，除了250 hPa溫度場(chǎng)的均方根誤差有所變壞，其他預(yù)報(bào)變量的均方根誤差都有所減少。東亞地區(qū)各預(yù)報(bào)量的距平相關(guān)系數(shù)與控制試驗(yàn)相當(dāng)，溫度場(chǎng)和高度場(chǎng)的均方根誤差有所減少。南半球基本持平。

圖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月模式預(yù)報(bào)變量評(píng)分卡 (紅色：變好,綠色：變差,灰色：相當(dāng))Fig.7 Model prediction variable score card of July 2017 (red: better, green: worse, gray: equal)

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

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

5 結(jié)論與討論

GRAPES_GFS2.0的垂直速度在模式高層出現(xiàn)計(jì)算噪音會(huì)導(dǎo)致積分不穩(wěn)定，甚至積分中斷。去掉模式物理過程和模式地形后這個(gè)噪音仍然存在，因此認(rèn)為動(dòng)力框架計(jì)算誤差是計(jì)算噪音產(chǎn)生的主要原因。

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