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

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

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

1 引言

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

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

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

為了分析解決半隱式方案中t+Δt時刻非線性項顯式處理帶來的計算不穩(wěn)定，Cullen（2001）將預(yù)估-校正（Predictor-Corrector，簡稱：P-C）方法用于二時間層淺水波模式來提高模式的積分穩(wěn)定性。預(yù)估-校正法一般將原來一步完成的積分分成多步計算，多步計算可以是全顯式計算，比如Kurihara（1965）在預(yù)估（P 階段）采用蛙躍格式，在校正（C 階段）采用二階精度的時間平均法；可以是半隱式處理，比如Kar（2012）在預(yù)估-校正的每個階段都將一個三階的Runge-Kutta 時間差分方案與一個重力波項二階精度時間平均法進(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ù)，降低模式的計算精度或減少積分時間步長，降低模式計算效率。本研究將在CMA-GFS 上發(fā)展預(yù)估-校正的半隱式半拉格朗日方案（SISL/P-C），以減少時間外插的影響提高模式的計算精度和穩(wěn)定性。

2 SISL/P-C 方案

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

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

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

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

其他預(yù)報變量則可以寫成關(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）右端項的函數(shù)，是時間變化項。上述方程組在求解過程中是作為已知項放在方程組的右端顯式表達(dá)，在傳統(tǒng)的二時間層SISL 方案中是通過時間外插得到的

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

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

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

預(yù)估P：

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

3 理想試驗

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

3.1 Rossby-Haurhawtz 波試驗

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

圖1 Rossby-Haurhawtz 波試驗積分15 d 后700 hPa 高度場（a.對照試驗，b.SISL/P-C 試驗；單位：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 波試驗 （40°S—40°N，10°—110°E）積分15 d 后700 hPa 高度場 （黑色：初始值， 紅色：對照試驗，綠色：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)生羅斯貝波，本試驗就是模擬羅斯貝波的產(chǎn)生、演變和消亡過程。具體的試驗設(shè)置參照J(rèn)ablonowski 等（2008）第5 個干動力框架標(biāo)準(zhǔn)試驗。模式水平分辨率0.5°，垂直均勻分成60 層，每層厚400 m，時間積分步長分別為600 s 和1200 s。圖3 是模式積分15 d 的700 hPa 高度場。2 個試驗在時間步長分別取600 s 和1200 s 時都可以很清晰地模擬出山脈激發(fā)的羅斯貝波的激發(fā)和傳播以及消亡過程。氣流過山后在背風(fēng)坡都形成了明顯的地形強(qiáng)迫羅斯貝波。

圖3 地形激發(fā)的羅斯貝波傳播15 d 后700 hPa 高度場（a1、a2：對照試驗，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）

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

圖4 地形激發(fā)的羅斯貝波傳播15 d 后700 hPa（EQ—70°N，80°—160°E）高度場（實線：600 s，虛線：1200 s；黑色：對照試驗，紅色：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 平均地面氣壓變化（黑色：對照試驗，紅色：SISL/P-C；實線：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 斜壓波試驗

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

圖6 斜壓波試驗?zāi)Ｊ椒e分9 天850 hPa 上的溫度場（單位：K）（a.對照試驗，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 斜壓波試驗?zāi)Ｊ椒e分9 天850 hPa（30—70°N，120°E —120°W）溫度場（a.對照試驗，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 實際資料試驗

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

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

圖8 2018 年7 月12 日12 時（世界時）起報72 h 預(yù)報的500 hPa 高度場（黑色線：SISL/P-C，紅色線：對照試驗；實線：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）

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

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

圖10 SISL/P-C 預(yù)報試驗的綜合評分卡（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 平均的全球平均海平面氣壓隨時間的變化。圖中黑色實線是對照試驗的結(jié)果，8 d 預(yù)報的平均地面氣壓下降0.65 hPa。SISL/P-C 8 d 預(yù)報的平均地面氣壓與初始平均地面氣壓相比損失都有所減少，當(dāng)時間步長與對照試驗一樣為300 s 時，8 d 預(yù)報平均地面氣壓下降0.5 hPa，比對照試驗明顯減少。

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

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

5 結(jié)論與討論

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

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

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

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

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