預(yù)估-校正的半隱式半拉格朗日時間積分方案及其在CMA-GFS 模式中的應(yīng)用*

張紅亮 沈?qū)W順 蘇 勇

中國氣象局地球系統(tǒng)數(shù)值預(yù)報(bào)中心，北京，100081

中國氣象科學(xué)研究院災(zāi)害天氣國家重點(diǎn)實(shí)驗(yàn)室，北京，100081

1 引言

半隱式-半拉格朗日（SISL）方法是Robert（1981）提出，并在1985 年用于靜力原始方程模式中（Robert，et al，1985）。由于SISL 方法不受Courant-Friedrichs-Lewy（CFL）的限制，在大時間步長下仍然可以保證模式的穩(wěn)定和精確，很快就被應(yīng)用到了業(yè)務(wù)數(shù)值天氣預(yù)報(bào)模式中。歐洲中期天氣預(yù)報(bào)中心（ECMWF）在1991 年發(fā)展了基于三時間層SISL 方法的業(yè)務(wù)靜力譜模式（Ritchie,et al,1995），1996 年又用二時間層SISL 模式替換了三時間層SISL 模式，計(jì)算效率提高了1 倍（Temperton，et al，2001）。此外，英國氣象局2005 年（Davies，et al，2005）、加拿大氣象 局2002 年（Yeh，et al，2002）、中國氣象局2009 年（薛紀(jì)善等，2008；陳德輝等，2006；沈?qū)W順等，2017）都開發(fā)了基于二時間層SISL 的全球非靜力格點(diǎn)模式。

二時間層SISL 模式假設(shè)模式格點(diǎn)上的空氣質(zhì)點(diǎn)在t到t+Δt的一個時間步長內(nèi)沿著拉格朗日軌跡做勻速直線運(yùn)動，移動速度就是軌跡中間點(diǎn)時刻的速度，用t和t?Δt時刻速度線性外插得到。隨著模式分辨率的提高，在急流軸等風(fēng)速梯度大值區(qū)時間線性外插會造成拉格朗日軌跡嚴(yán)重的計(jì)算誤差，產(chǎn)生計(jì)算噪音，被認(rèn)為是造成二時間層SISL 模式積分不穩(wěn)定的主要原因。Hortal（2002）將勻速直線運(yùn)動假定變成勻加速直線運(yùn)動假定，將運(yùn)動方程在拉格朗日上游點(diǎn)展開，保留二階導(dǎo)數(shù)項(xiàng)，構(gòu)造了穩(wěn)定外插的二時間層方案（Stable extrapolation two time level scheme,SETTLS），減少了急流軸附近的虛假計(jì)算噪音。SETTLS 方案本質(zhì)上仍然是時間外插方案，不能從根本上解決時間外插帶來的計(jì)算不穩(wěn)定問題（Durran，2004）。

二時間層SISL 模式通過引入?yún)⒖即髿鈱⒎匠逃叶藲鈮禾荻攘?、科里奧利力和浮力項(xiàng)分成代表快波的線性項(xiàng)（L）和代表慢波的非線性項(xiàng)（N），因?yàn)橄嗨俣鹊牟煌?，t+Δt時刻的L項(xiàng)采取隱式處理以增強(qiáng)穩(wěn)定性，其他項(xiàng)采用顯式處理。t+Δt時刻的N項(xiàng)也是顯式表達(dá)，用t和t?Δt時刻N(yùn)項(xiàng)線性外插得到，但在實(shí)踐過程中這種線性外插會造成地形附近的計(jì)算噪音（Coiffier，et al，1987）。

為了分析解決半隱式方案中t+Δt時刻非線性項(xiàng)顯式處理帶來的計(jì)算不穩(wěn)定，Cullen（2001）將預(yù)估-校正（Predictor-Corrector，簡稱：P-C）方法用于二時間層淺水波模式來提高模式的積分穩(wěn)定性。預(yù)估-校正法一般將原來一步完成的積分分成多步計(jì)算，多步計(jì)算可以是全顯式計(jì)算，比如Kurihara（1965）在預(yù)估（P 階段）采用蛙躍格式，在校正（C 階段）采用二階精度的時間平均法；可以是半隱式處理，比如Kar（2012）在預(yù)估-校正的每個階段都將一個三階的Runge-Kutta 時間差分方案與一個重力波項(xiàng)二階精度時間平均法進(jìn)行組合。最早在二時間層SISL 業(yè)務(wù)模式中使用預(yù)估-校正時間積分方案的是加拿大氣象局的GEM 模式（C?te，et al，1998）。后來Cullen（2001）將此方案用于ECMWF 的IFS模 式 中，Diamantakis 等（2007）在UKMO 的New Dynamical Core 上，Wood 等（2014）在ENDGAME中都使用了這個時間積分方案。

CMA-GFS 采用的是二時間層SISL 時間積分方案同樣面臨時間外插造成的模式積分不穩(wěn)定問題。通常解決的辦法是增大半隱式系數(shù)，降低模式的計(jì)算精度或減少積分時間步長，降低模式計(jì)算效率。本研究將在CMA-GFS 上發(fā)展預(yù)估-校正的半隱式半拉格朗日方案（SISL/P-C），以減少時間外插的影響提高模式的計(jì)算精度和穩(wěn)定性。

2 SISL/P-C 方案

CMA-GFS 采用的二時間層SISL 通用表達(dá)式如下

式中，X為三維速度場V、擾動無量綱氣壓П′、擾動位溫 θ′等預(yù)報(bào)變量；L為引入?yún)⒖即髿庵蠓匠探M右端的線性項(xiàng)；N為方程組右端扣除線性項(xiàng)之后的殘差小項(xiàng)，也就是非線性項(xiàng)；α是半隱式權(quán)重系數(shù)；A是拉格朗日軌跡到達(dá)點(diǎn)，也是模式格點(diǎn)；D是到達(dá)A點(diǎn)的拉格朗日上游點(diǎn)。

將式（1）重新寫成線性方程形式

方程組（2）通過消元法得到關(guān)于t+?t時刻П′的Helmholtz 方程組

其他預(yù)報(bào)變量則可以寫成關(guān)于П′的顯式表達(dá)式

方程組的推導(dǎo)和系數(shù)表達(dá)式參見薛紀(jì)善等（2008）。系數(shù)ξu1，ξu2，ξu3，ξv1，ξv2，ξv3···等是關(guān)于?t、α、dφ等的常數(shù)，ξu0，ξv0，ξw0，AΠ···是式（2）右端項(xiàng)的函數(shù)，是時間變化項(xiàng)。上述方程組在求解過程中是作為已知項(xiàng)放在方程組的右端顯式表達(dá)，在傳統(tǒng)的二時間層SISL 方案中是通過時間外插得到的

同樣，半拉格朗日軌跡上的質(zhì)點(diǎn)移動速度也是采用時間外插得到

時間外插在梯度大的地方會造成計(jì)算噪音，導(dǎo)致模式積分不穩(wěn)定。

為了減少時間外插的影響，CMA-GFS 在現(xiàn)有的傳統(tǒng)二時間層SISL 方案的基礎(chǔ)上開發(fā)了SISL/P-C 方案，式（2）通過兩步來實(shí)現(xiàn)。

預(yù)估P：

式（10）、（11）表明CMA-GFS 動力過程需要在預(yù)估P 和校正C 計(jì)算兩次，也就是方程（3）—（7）要求解兩次。在預(yù)估P 過程中拉格朗日軌跡中間點(diǎn)的速度可以用式（9）計(jì)算或者直接用到達(dá)點(diǎn)t時刻的速度代替（穩(wěn)定性更好），估算出拉格朗日上游點(diǎn)的位置D；由式（8）得到。預(yù)估P 完成之后預(yù)報(bào)變量有了初猜場式（8）和（9）將變成式（12）和（13）

3 理想試驗(yàn)

DCMIP（Dynamical Core Model Intercomparison Project）（Jablonowski，et al，2008）是一個通過設(shè)計(jì)一系列理想試驗(yàn)比較驗(yàn)證模式動力框架性能的計(jì)劃。文中用DCMIP 中的3 個經(jīng)典理想試驗(yàn)：Rossby-Haurhawtz 波、地形激發(fā)的羅斯貝波和斜壓波驗(yàn)證SISL/P-C 方案的精確性、守恒性和穩(wěn)定性。

3.1 Rossby-Haurhawtz 波試驗(yàn)

Rossby-Haurhawtz 波是正壓渦度方程的解析解，用來模擬天氣尺度波動的傳播和演變過程（Thuburn，et al，2000）。假設(shè)在赤道上放置4 個自西向東傳播的波，理想情況下，這4個波在傳播過程中會一直保持形狀不發(fā)生衰減。這個試驗(yàn)可以展示動力框架在積分過程中的耗散、質(zhì)量守恒和相速度誤差等性質(zhì)。具體試驗(yàn)設(shè)置參考Jablonowski等（2008）第4 個標(biāo)準(zhǔn)試驗(yàn)。模式水平分辨率1.0°，垂直不等距分為60 層，時間步長1200s，模式積分15 d 的700 hPa 高度場如圖1 所示。對照試驗(yàn)和SISL/P-C 試驗(yàn)都可以模擬出4個波自西向東的傳播過程，位相傳播速度基本一致。在圖1 中截取一個完整的波（40°S—40°N，10°—110°E），并與初始場中一個完整波（40°S—40°N，38.25°—138.25°E）進(jìn)行疊加如圖2 所示。與初始場相比，對照試驗(yàn)積分15 d 后的波形有了很大的衰減，采用SISL/P-C方案后波形與初值基本上保持一致，只是南北方向略有衰減。

圖1 Rossby-Haurhawtz 波試驗(yàn)積分15 d 后700 hPa 高度場（a.對照試驗(yàn)，b.SISL/P-C 試驗(yàn)；單位：gpm）Fig.1 700 hPa height in the Rossby-Haurhawtz wave experiment after 15 d of integration（a.control，b.SISL/P-C；unit：gpm）

圖2 Rossby-Haurhawtz 波試驗(yàn) （40°S—40°N，10°—110°E）積分15 d 后700 hPa 高度場 （黑色：初始值， 紅色：對照試驗(yàn)，綠色：SISL/P-C ）Fig.2 700 hPa height in the Rossby-Haurhawtz wave experiment（40°S—40°N，10°—110°E）after 15 d ofintegration（black：initial value，red：control，green：SISL/P-C）

3.2 地形激發(fā)的羅斯貝波

假設(shè)在北半球中緯度地區(qū)放置一個鐘形地形，流經(jīng)地形的平直西風(fēng)氣流在地形的激發(fā)下會產(chǎn)生羅斯貝波，本試驗(yàn)就是模擬羅斯貝波的產(chǎn)生、演變和消亡過程。具體的試驗(yàn)設(shè)置參照J(rèn)ablonowski 等（2008）第5 個干動力框架標(biāo)準(zhǔn)試驗(yàn)。模式水平分辨率0.5°，垂直均勻分成60 層，每層厚400 m，時間積分步長分別為600 s 和1200 s。圖3 是模式積分15 d 的700 hPa 高度場。2 個試驗(yàn)在時間步長分別取600 s 和1200 s 時都可以很清晰地模擬出山脈激發(fā)的羅斯貝波的激發(fā)和傳播以及消亡過程。氣流過山后在背風(fēng)坡都形成了明顯的地形強(qiáng)迫羅斯貝波。

圖3 地形激發(fā)的羅斯貝波傳播15 d 后700 hPa 高度場（a1、a2：對照試驗(yàn)，b1、b2.SISL/P-C；a1、b1.1200 s，a2、b2.600 s；單位：gpm）Fig.3 700 hPa height in the mountain-induced Rossby wave experiment after 15 d of integration（a1，a2.control，b1，b2.SISL/P-C；a1，b1.1200 s，a2，b2.600 s；unit：gpm）

對照試驗(yàn)中模擬的山脈下游的槽、脊強(qiáng)度與其他模式相比較弱（蘇勇等，2018）。本試驗(yàn)中選取該區(qū)域進(jìn)一步仔細(xì)比較兩組試驗(yàn)的差別，如圖4 所示。不論時間步長取600 s 或是1200 s，SISL/P-C模擬的槽和脊強(qiáng)度都比對照試驗(yàn)加強(qiáng)，明顯緩解了系統(tǒng)偏弱的問題。圖5 給出了0—15 d 平均地面氣壓的變化。SISL/P-C 在2 組試驗(yàn)中經(jīng)過15 d 的預(yù)報(bào)，地面氣壓的變化明顯低于對照試驗(yàn)，當(dāng)時間步長取600 s 時，積分15 d 后對照試驗(yàn)中地面氣壓降低了 1.13 hPa，SISL/P-C 降低0.57 hPa，是對照試驗(yàn)的一半；當(dāng)時間步長取1200 s，積分15 d 后SISL/P-C的地面氣壓變化比對照試驗(yàn)減少0.2 hPa。

圖4 地形激發(fā)的羅斯貝波傳播15 d 后700 hPa（EQ—70°N，80°—160°E）高度場（實(shí)線：600 s，虛線：1200 s；黑色：對照試驗(yàn)，紅色：SISL/P-C；單位：gpm）Fig.4 700 hPa height（EQ—70°N，80° —160°E）in the mountain-induced Rossby wave experiment after 15 d of integration（solid：600 s，dashed：1200 s；black：control，red：SISL/P-C；unit：gpm）

圖5 地形激發(fā)的羅斯貝波傳播0—15 d 平均地面氣壓變化（黑色：對照試驗(yàn)，紅色：SISL/P-C；實(shí)線：600 s，虛線：1200 s；單位：hPa）Fig.5 Averaged surface pressure changes over 0—15 d in the mountain-induced Rossby wave experiment（black：control，red：SISL/P-C；solid line：600 s，dashed line：1200 s；unit：hPa）

3.3 斜壓波試驗(yàn)

斜壓波試驗(yàn)是在一個平直的西風(fēng)氣流上疊加一個擾動的緯向風(fēng)場后激發(fā)出了一個短時斜壓波，這個斜壓波在第9 天發(fā)展最為強(qiáng)盛，然后逐漸衰亡。具體擾動設(shè)置參見Jablonowski 等（2008）第1 個干動力框架標(biāo)準(zhǔn)試驗(yàn)。模式水平分辨率0.5°，垂直采用業(yè)務(wù)60 層不均勻分層，時間步長600 s。兩組試驗(yàn)都模擬出擾動風(fēng)場的發(fā)生、發(fā)展和加強(qiáng)過程，如圖6 所示。與Jablonowski 中T106 結(jié)果相比，積分第9 天兩組試驗(yàn)都激發(fā)出了3 個波峰，波峰位置和強(qiáng)度與T106 基本一致。進(jìn)一步將時間步長延長到1200 s 和1800 s，選擇子區(qū)域（30°—70°N，120°E —100°W）對比3 種時間步長下第9 天850 hPa溫度場，如圖7 所示。對照試驗(yàn)在600 s 積分時間步長下波的發(fā)展最強(qiáng)烈，積分步長增加1 倍波的振幅已經(jīng)有所衰減；增加2 倍時衰減已經(jīng)非常明顯。SISL/P-C 試驗(yàn)中3 種時間步長下第9 天的溫度場還是基本上重合在一起，沒有因?yàn)闀r間步長的增大而明顯衰減，說明SISL/P-C 方案有非常好的穩(wěn)定性。這也使得CMA-GFS 延長積分時間步長，提高計(jì)算效率成為可能。

圖6 斜壓波試驗(yàn)?zāi)Ｊ椒e分9 天850 hPa 上的溫度場（單位：K）（a.對照試驗(yàn)，b.SISL/P-C；單位：K）Fig.6 850 hPa temperature in the baroclinic wave experiment after 9 d of integration（a.control，b.SISL/P-C；unit：K）

圖7 斜壓波試驗(yàn)?zāi)Ｊ椒e分9 天850 hPa（30—70°N，120°E —120°W）溫度場（a.對照試驗(yàn)，b.SISL/P-C；黑色：600 s，藍(lán)色：1200 s，紅色：1800 s；單位：K）Fig.7 850 hPa temperature（30°—70°N，120°E—120°W）in the baroclinic wave experiment after 9 d of integration（a.control，b.SISL/ P-C；black：600 s，blue：1200 s，red：1800 s；unit：K）

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

CMA-GFS_V3.0（蘇勇等，2020）水平分辨率0.25°，垂直62 層，模式層頂位于36 km。物理過程包括長短波輻射、積云對流參數(shù)化過程、微物理降水過程、邊界層和陸面過程以及地形重力波等過程（沈?qū)W順等，2017；馬占山等，2016；宮宇等，2018）。采用NCEP 的FNL（Final operational global analysis data）資料作為模式的初始場，試驗(yàn)時段為2018 年7 月1 日到31 日，每天12 時（世界時）開始進(jìn)行8 d預(yù)報(bào)。對照試驗(yàn)中半隱式系數(shù)為0.72，SISL/P-C中是0.55，時間積分步長分別取300 s 和450 s。

圖8 是2018 年7 月12 日72 h 預(yù)報(bào)的500 hPa高度場。圖中黑色線是SISL/P-C 試驗(yàn)，紅色線是對照試驗(yàn)；實(shí)線和虛線代表時間步長分別為300 s和450 s。SISL/P-C 試驗(yàn)中兩種時間步長下預(yù)報(bào)的形勢場并沒有明顯差別，而且與對照試驗(yàn)在300 s時間步長下的預(yù)報(bào)結(jié)果基本吻合。當(dāng)對照試驗(yàn)中時間步長延長到450 s 時72 h 預(yù)報(bào)的形勢場有明顯衰減，這個結(jié)論與理想場試驗(yàn)的結(jié)果是一致的：SISL/P-C 方案是穩(wěn)定的，對照試驗(yàn)的穩(wěn)定性差。所以在0.25°的水平分辨率下對照試驗(yàn)的時間步長為300 s，SISL/P-C 試驗(yàn)的時間步長在此基礎(chǔ)上增加到450 s。

圖8 2018 年7 月12 日12 時（世界時）起報(bào)72 h 預(yù)報(bào)的500 hPa 高度場（黑色線：SISL/P-C，紅色線：對照試驗(yàn)；實(shí)線：300 s，虛線：450 s；單位：gpm）Fig.8 500 hPa height after 72 h of integration starting from 12：00 UTC 12 July 2018（solid：300 s，dash：450 s；black：SISL/P-C，red：control；unit：gpm）

兩組試驗(yàn)月平均的5 d 預(yù)報(bào)500 hPa 高度場的距平相關(guān)系數(shù)和預(yù)報(bào)偏差如圖9 所示。SISL/PC 試驗(yàn)中在兩組時間步長下500 hPa 高度場的距平相關(guān)系數(shù)和預(yù)報(bào)偏差沒有明顯區(qū)別。與對照試驗(yàn)相比，500 hPa 距平相關(guān)系數(shù)略有提高；預(yù)報(bào)偏差有明顯改進(jìn)。溫度場的檢驗(yàn)結(jié)果相同，風(fēng)場的檢驗(yàn)基本持平，沒有明顯改進(jìn)。

圖10 是SISL/P-C 預(yù)報(bào)試驗(yàn)的綜合評分卡，分東亞、北半球、南半球、赤道4 個區(qū)域分別對風(fēng)場、溫度場、高度場進(jìn)行檢驗(yàn)，左列為預(yù)報(bào)時效1—8 d的距平相關(guān)系數(shù)，右列為均方根誤差，紅色表示正效果，綠色表示負(fù)效果，灰色為中性，箭頭越大效果越顯著。SISL/P-C 的時間步長與對照試驗(yàn)一樣為300 s 時，模式的整體預(yù)報(bào)性能都有所改善，特別是均方根誤差有明顯降低，南半球的各個預(yù)報(bào)場以及熱帶地區(qū)的風(fēng)場的距平相關(guān)系數(shù)也有明顯提高。SISL/P-C 的時間步長提高到450 s，除了東亞和北半球地區(qū)溫度場的距平相關(guān)系數(shù)有負(fù)影響，其他地區(qū)是中性偏正的影響。

圖10 SISL/P-C 預(yù)報(bào)試驗(yàn)的綜合評分卡（a.300 s，b.450 s）Fig.10 Comprehensive scoring cards for prediction experiment of SISL/P-C（a.300 s，b.450 s）

圖11 是2018 年7 月31 d 平均的全球平均海平面氣壓隨時間的變化。圖中黑色實(shí)線是對照試驗(yàn)的結(jié)果，8 d 預(yù)報(bào)的平均地面氣壓下降0.65 hPa。SISL/P-C 8 d 預(yù)報(bào)的平均地面氣壓與初始平均地面氣壓相比損失都有所減少，當(dāng)時間步長與對照試驗(yàn)一樣為300 s 時，8 d 預(yù)報(bào)平均地面氣壓下降0.5 hPa，比對照試驗(yàn)明顯減少。

圖11 2018 年7 月全球平均的海平面氣壓隨時間的變化Fig.11 Changes of global mean sea level pressure over time in July 2018

SISL/P-C 的預(yù)估、校正過程中動力框架要積分兩次，第2 次積分過程中t時刻格點(diǎn)上的值不再重復(fù)計(jì)算，并且將預(yù)估值作為求解Helmholtz 方程組的初猜場極大減少了方程組的余差和迭代收斂次數(shù)，在與對照試驗(yàn)相同的時間步長下8 d 預(yù)報(bào)時間僅增加了40%。由于SISL/P-C 方案有非常好的計(jì)算穩(wěn)定性，0.25°分辨率下時間步長可以延長50%到450 s，模式預(yù)報(bào)結(jié)果和效果仍然是穩(wěn)定的，8 d 預(yù)報(bào)時間減少20%，計(jì)算效率明顯提高。

5 結(jié)論與討論

為了減少二時間層SISL 方案中時間外插對預(yù)報(bào)穩(wěn)定性和精度的影響，在CMA-GFS_V3.0 的基礎(chǔ)上發(fā)展了SISL/P-C 方案。SISL/P-C 中動力框架要完整地計(jì)算兩次，為了最大限度地保證預(yù)估P 過程的精度，預(yù)估P 過程與CMA-GFS_V3.0 一致；在校正C 過程中原來需要時間外插的拉格朗日上游點(diǎn)移動速度根據(jù)預(yù)估值可以線性內(nèi)插得到，需要顯式處理的t+?t時刻非線性項(xiàng)也可以根據(jù)預(yù)估值直接計(jì)算。

文中用3 個經(jīng)典理想試驗(yàn)（Rossby-Haurhawtz波、地形激發(fā)的羅斯貝波和斜壓波試驗(yàn)）證明了相比于對照試驗(yàn)，SISL/P-C 的波形衰減明顯減少，計(jì)算精度顯著提高。將積分時間步長延長2—3 倍后SISL/P-C 長時間積分后波的形狀、強(qiáng)度和移動速度仍然保持不變，模式的穩(wěn)定性得到了很大的提高，這也為實(shí)際預(yù)報(bào)中增大積分步長提高計(jì)算效率提供了支持和保障。

本研究也進(jìn)行了1 個月的實(shí)際資料的對比試驗(yàn)，SISL/P-C 的總體預(yù)報(bào)效果優(yōu)于對照試驗(yàn)。不同積分步長下SISL/P-C 72 h 預(yù)報(bào)的500 hPa 高度場和溫度場基本上是一樣的，對照試驗(yàn)隨著時間步長的增大有明顯衰減，新方案提高了模式的積分穩(wěn)定性。

基于SISL/P-C 方案等改進(jìn)生成的CMA-GFS_V3.1 四維變分同化系統(tǒng)2020 年汛期前正式業(yè)務(wù)運(yùn)行。模式水平分辨率0.25 °,垂直不等距分為87 層，模式層頂延伸到0.1 hPa，時間積分步長由300 s 提高到450 s。

目前，開發(fā)的SISL/P-C 方案是CMA-GFS_V3.1動力框架的計(jì)算方案。Wood 等（2014）和陳子通等（2020）在預(yù)估-校正方案中增加了物理過程的影響。采用SISL/P-C 方案后CMA-GFS 動力框架的半隱式系數(shù)由0.72 減少到0.55，動力框架近于二階計(jì)算精度，但是物理過程和動力框架的耦合還是一階精度，所以進(jìn)一步提高CMA-GFS_V3.1 的動力-物理耦合精度將是未來的工作重點(diǎn)。