引言

二时间层半隐式半拉格朗日时间积分方案(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中将当前时刻到达点的速度作为质点的移动速度，避免了时间外插带来的计算噪音，增加了模式的计算稳定性。

另一方面，拉格朗日上游点一般不在模式格点上，上游点上的物理量可以通过周围网格点高阶插值得到，比如Hermite三次插值方案、cubic拉格朗日插值方案和准三次cubic拉格朗日插值方案，以及水平三阶、垂直五阶的插值方案等。随着模式分辨率越来越高, 插值误差对模式预报误差的影响越来越明显(Polichtchouk et al, 2020)。Ritchie(1986)早就提出了在拉格朗日平流的基础上通过减少插值来减少平滑对计算精度的影响。Rańič (1995)提出非插值的半拉格朗日平流方案可以减少模式的耗散误差。英国气象局的New Dynamical Core和ECMWF的IFS模式都在位温方程中引入垂直非插值拉格朗日平流方案，减少了垂直插值对位温造成的计算误差，提高了模式的计算精度(Diamantakis et al, 2007; Ritchie et al, 1995)。

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展开，得到：

$ \frac{{X_A^{t + \Delta t} - X_D^t}}{{\Delta t}} = \alpha \left({L_A^{t + \Delta t} + N_A^{t + \Delta t}} \right) + (1 - \alpha)\left({L_D^t + N_D^t} \right) $

(1)

式中：A表示t+Δt时刻拉格朗日到达点；D表示t时刻拉格朗日上游点；L是引入参考大气之后方程右端的线性项，N是扣除线性项的残差项，也是非线性项，L和N的具体表达式参见薛纪善和陈德辉(2008)、沈学顺等(2017)；α是半隐式时间系数，为了保持积分稳定性，GRAPES_GFS中1>α>0.5；N_A^t+Δt是t+Δt到达点的非线性项，通过t和t-Δt时刻的非线性项时间线性外插得到：

$ N_A^{t + \Delta t} = 2N_A^t - N_A^{t - \Delta t} $

(2)

线性方程组(1)的求解最后转为求解关于无量纲扰动气压的Helmholtz方程组。拉格朗日轨迹上游点计算在极区采用Mcdonald and Bates(1989)方法，在中低纬度采用的是Ritchie(1987；1988)方法。拉格朗日轨迹中间点速度为：

$ {\mathit{\boldsymbol{V}}^{t + \frac{1}{2}\Delta t}} = \frac{3}{2}{\mathit{\boldsymbol{V}}^t} - \frac{1}{2}{\mathit{\boldsymbol{V}}^{t - \Delta t}} $

(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.2 拉格朗日上游点计算的影响

本节通过敏感性试验，讨论拉格朗日上游点计算对计算噪音的影响。

从拉格朗日到达点出发，沿着i, j, k方向分别计算一个时间步长下质点反向移动的距离得到拉格朗日上游点的位置。GRAPES_GFS的4套网格点计算4套拉格朗日轨迹并在上游点上进行插值。4套上游点的计算方案和插值方案是一样的，如果是拉格朗日轨迹计算导致的上游点位置计算误差带来的计算噪音，则4套上游点的位置计算都会产生计算噪音。但是上游点上不同变量的插值误差对模式预报精度的影响是不一样的。

试验第一步：讨论拉格朗日上游点位置计算误差对计算噪音的敏感性。分别设置4套拉格朗日平流在i, j, k方向的移动速度为零，在此方向上上游点和到达点处在同一个网格点上，在此方向上也不需要进行空间插值。12组试验中，只有将垂直速度和位温所在格点的垂直拉格朗日平流速度设为零时，可以消除计算噪音(图 4)，其他设置下，计算噪音仍然存在。试验结果表明计算噪音的产生与拉格朗日上游点位置计算无关，垂直速度或位温在垂直方向的插值误差是造成计算噪音的主要原因。

试验第二步：拉格朗日轨迹计算方案不变，垂直速度和位温在拉格朗日上游点上分别只进行水平插值，不进行垂直插值，以讨论垂直速度和位温的垂直插值误差对计算噪音的影响。试验结果表明拉格朗日上游点上的位温，如果不进行垂直插值而只进行水平插值，计算噪音也会消除(图 5a)；但是，垂直速度不管是否进行垂直插值，计算噪音都会存在。图 5b是垂直速度在上游点上不进行垂直插值的情况，可以看到计算噪音仍然存在，最大上升速度为7 m·s^-1。

试验第三步：在第二步试验的基础上，将位温的拉格朗日上游点在垂直方向上近似到最接近的模式层，拉格朗日上游点上只需要进行水平插值，不做垂直插值。试验结果表明，垂直速度场是光滑的，计算噪音没有出现(图 5c)。

以上三步试验表明，拉格朗日上游点上位温的垂直插值会造成模式高层垂直速度的计算噪音。

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

第2节分析认为，位势温度及其相关项在拉格朗日上游点上的垂直插值是模式计算噪音产生的原因。本文据此对位温方程的拉格朗日垂直平流进行改造, 将垂直方向的拉格朗日上游点近似到最接近的模式面, 垂直方向不再进行插值。

在没有源汇项的情况下，位势温度θ的热量方程如下：

$ \frac{{{\rm{d}}\theta }}{{{\rm{d}}t}} = 0 $

(4)

引入满足静力平衡的一维参考大气廓线θ(z)，得到扰动位势温度θ′=θ-θ(z)，式(4)变成：

$ \frac{{{\rm{d}}\theta }}{{{\rm{d}}t}} = \frac{{{\rm{d}}{\theta ^\prime }}}{{{\rm{d}}t}} + w\frac{{\partial \bar \theta }}{{\partial z}} = 0 $

(5)

将$\frac{{{\rm{d}}{\theta ^\prime }}}{{{\rm{d}}t}}$进行欧拉展开得到：

$ \frac{{{\rm{d}}{\theta ^\prime }}}{{{\rm{d}}t}} = \frac{{\partial {\theta ^\prime }}}{{\partial t}} + u\frac{{\partial {\theta ^\prime }}}{{a\cos \varphi \partial \lambda }} + v\frac{{\partial {\theta ^\prime }}}{{a\partial \varphi }} + \hat w\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}} + w\frac{{\partial \bar \theta }}{{\partial z}} = 0 $

(6)

式中：u, v, ŵ分别是模式面水平和垂直风速，w是高度面上的垂直风速。ẑ是模式面高度，z是海拔高度。

对该方程中的垂直平流项进行修正(Diamantakis et al, 2007)，假设r_A是θ在t+Δt时刻所在拉格朗日到达点，在一个时间步长内它沿着拉格朗日轨迹以(u, v, ŵ)做反向匀速直线运动到达拉格朗日上游点r_D。为了避免位温及其相关项在垂直方向的插值，可以将r_D的垂直分量近似到最近的模式层r_Dl，近似的垂直速度ŵ^*等见式(7)

$ \begin{array}{l} {w^*} = \frac{{{r_A} - {r_{Dl}}}}{{{\rm{d}}t}}\\ \hat w = \frac{{{r_A} - {r_D}}}{{{\rm{d}}t}}\\ {r_{Dl}} = n{\mathop{\rm int}} \left({{r_D}} \right) \end{array} $

(7)

式中nint是取整函数。

位温方程可以重新写成

$ \frac{{{\theta ^{\prime t + \Delta t}} - \theta _{Dl}^{\prime t}}}{{\Delta t}} = - \left({\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}} - w\frac{{\partial \bar \theta }}{{\partial z}} $

(8)

对比式(8)和式(6)，不难看出公式修改后只是在方程右端增加了一个非线性项$ - \left({\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}}$。对式(8)进行2TSISL展开如下：

$ {\begin{array}{*{20}{c}} {{\theta ^\prime }_A^{t + \Delta t} = {\theta ^\prime }_{Dl}^t - \Delta t\left[ {\left( {\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}} + w\frac{{\partial \bar \theta }}{{\partial z}}} \right]_{Dl}^t - }\\ {\Delta t\left[ {\left( {\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}} + w\frac{{\partial \bar \theta }}{{\partial z}}} \right]_A^{t + \Delta t}} \end{array}} $

(9)

即：

$ {\theta _A^{\prime t + \Delta t} = - \Delta t\left({w\frac{{\partial \bar \theta }}{{\partial z}}} \right)_A^{t + \Delta t} + }\\ {\begin{array}{*{20}{c}} {\left\{ {{\theta ^\prime } - \Delta t\left[ {\left({\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}} + w\frac{{\partial \bar \theta }}{{\partial z}}} \right]} \right\}_{Dl}^t - }\\ {\left[ {\Delta t\left({\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}}} \right]_A^{t + \Delta t}} \end{array}} $

(10)

其中${\left[ {\Delta t\left({\hat w - {{\hat w}^*}} \right)\frac{{\partial {\theta ^\prime }}}{{\partial \hat z}}} \right]_A^{t + \Delta t}}$可由t和t-Δt的值通过线性外插得到。

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温度场的均方根误差有所变坏，其他预报变量的均方根误差都有所减少。东亚地区各预报量的距平相关系数与控制试验相当，温度场和高度场的均方根误差有所减少。南半球基本持平。

图 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的变化量。

5 结论与讨论

GRAPES_GFS2.0的垂直速度在模式高层出现计算噪音会导致积分不稳定，甚至积分中断。去掉模式物理过程和模式地形后这个噪音仍然存在，因此认为动力框架计算误差是计算噪音产生的主要原因。

GRAPES_GFS动力框架采用的是2TSISL方案。拉格朗日轨迹计算包括拉格朗日上游点的位置计算和上游点上的插值。诊断分析表明：拉格朗日上游点位置计算以及其他变量在拉格朗日上游点上的水平和垂直插值，都不会产生垂直速度场的计算噪音，只有位势温度在拉格朗日上游点处的垂直插值可以造成垂直速度的计算噪音。在位温方程中引入垂直非插值拉格朗日平流方案后，将位势温度的拉格朗日上游点在垂直方向上近似到最接近的模式层，从而避免了位温在垂直方向的插值带来的计算误差。2013年7月的实际资料试验表明，新的计算方案可以消除模式中高层的计算噪音，并且可以提高模式在北半球和热带地区的预报性能，减少预报的均方根误差。8 d预报质量损失严重的现象得到了明显的缓解。目前该方案已经在GRAPES_GFS2.0和GRAPES_GFS3.0上业务运行。