1 引 言

半隐式半拉格朗日(Semi-Implicit Semi-Lagrangian，以下简称SISL)时间积分方案具有绝对稳定性，即理论上时间步长的选取只取决于方案的计算精度，而不用考虑时间稳定性(Gravel et al，1993)，时间步长可以取得比显式的欧拉方案长得多，自从SISL时间积分方案提出以来在大气模式中得到了广泛应用。在应用SISL时间积分方案时，空气质点的轨迹计算对方案的精度和稳定性有重要影响，实际上轨迹计算的精度在很大程度上决定了模式计算的精度和稳定性(Staniforth et al，1991)。Ritche等(1994)假定空气质点在拉格朗日轨迹上是等速运动，采取中点近似，使轨迹中间点的速度近似等于半时间步时(t+Δt/2)的速度，再利用t和t-Δt时刻格点上的速度进行外插，得到网格点上t+Δt/2时刻的速度，然后进行空间插值迭代求解出上游点的位置。由于Ritche-Beaudoin上游点方案中引入了时间外插，因而会给SISL时间积分方案带来不稳定性。如果质点速度随轨迹变化，即不再是等速运动，当CFL(Coutant-Friedrichs-Lewy)计算稳定性判据大于1时，SISL方案就会不稳定(Durran，1999)。为了解决不稳定的问题，Hortal(2002)提出了两时间层稳定外插方案(Stable Extrapolation Two-Time-Level Scheme，以下简称SETTLS)，该方案假定质点沿拉格朗日轨迹做匀加速运动，这时轨迹时间中间点的位置不再是上游点与到达点的平均，从而改善了方案的稳定性。因为在SISL时间积分方案中对非线性项的处理也要涉及时间外插，因此，也可以将SETTLS方案用于处理非线性项，使得非线性项的计算在等权重时达到二阶精度。

在SISL时间积分方案中对线性项的处理也要慎重考虑。在半隐式时间积分方案中，线性项是对模式方程依照参考大气廓线进行线性化而得到的(Robert et al，1972)，一般对线性项进行非中央权重的时间平均，其能够有效地抑制虚假的半拉格朗日地形共振(Rivest et al，1994)，但同时也使得时间差分的精度由二阶变为一阶，而且，对大气真实的低波数波动振幅产生了较强的衰减(McDonald et al，1993; Bates et al，1993)。Simmons等(1978)也指出，线性化模式方程时参考大气廓线不能任意选择，否则可能造成时间积分不稳定。

中国气象局的全球和区域同化预报多尺度统一模式(GRAPES)采用了SISL时间积分方案，其轨迹计算使用Ritche等(1994)方案。模式积分有时会发生计算不稳定，模式同时又在时间积分方案上使用了较大的非中央权重系数处理相关项，这时模式时间积分方案只有一阶精度且对波动有显著的衰减。本研究拟将SETTLS时间积分方案引入GRAPES区域模式拉格朗日轨迹和非线性项的计算，同时对参考大气和线性项进行处理，将线性项的时间离散精度由一阶提高至二阶，并减小对波动的衰减作用。 2 上游点的SETTLS计算方案

按照Hortal(2002)提出的SETTLS，考虑到上游点在拉格朗日轨迹上可以通过求解拉格朗日轨迹方程得到。拉格朗日轨迹方程

式中，R和V 为三维的位置与速度矢量。

为了获得二阶精度的方案，将n+1时刻到达点的位置矢量 R ⁿ⁺¹_A在上游点附近做泰勒展开，保留二次项

式中，下标AV代表沿着拉格朗日轨迹上的某个平均，下标D和A分别表示上游点和到达点。将拉格朗日轨迹方程代入式(2)得上述方程代表了一个匀加速运动。

平均加速度可取为

将式(4)代入式(3)得可以证明此方案是线性稳定的(Durran，1999；Durran et al，2004)。该方案需要迭代求解，设取(dR)⁽⁰⁾=0，则有式中，k为迭代次数，亦可将此式表示为式中，实际上该方案描述了一个匀加速运动，拉格朗日轨迹的时间中间点不再是上游点和到达点的平均。上述方程是一个3维的矢量方程，需在球面上求解。

考虑球面上的算法，上游点可由以下迭代方案进行计算

式中，a为地球半径，r _D、r _A分别是上游点和到达点到地心的距离，φ为 r _D与 r _A的夹角；R^DA为上游点到到达点的旋转算子式中，φ、λ分别为纬度和经度，且p²+q²=1。

如果φ是小角度近似，则有

因此，于是就有在垂直方向对于第1次迭代，设上游点(D)为到达点(A)，φ₀=| V n_A|，下一步迭代，将(2 V(n)- V(n-1))和(2(n)-(n-1))插到第1步计算出的D(φ⁽¹⁾_D，λ⁽¹⁾_D，⁽¹⁾_D)点，如此循环，一般进行3次迭代就能够达到收敛精度。

可将 V (u，v)写为一般形式 V n_A+ [2 V_D^k(n)- V ^k_D(n-1)]。模式实际计算中是将矢量写成分量形式

与Ritche等(1994)方案相比，上述方案不需要首先计算中间点的位置，再计算得到上游点的位置，而是直接求得上游点的位置。 3 非线性项SETTLS计算方案与二阶精度非中央权重的线性项计算方案

目前GRAPES区域模式的SISL时间积分方案为

式中，权重系数ε取为0.7。该方案由于对线性项和非线性项都取非中央权重系数，只有一阶精度，而且会对波动产生较大的衰减作用，ε越接近于1，衰减作用越强。

如上所述，SISL时间积分方案中对线性项采取的是隐式处理，这样可以抑制快波项对时间步长的限制，使得SISL时间积分方案无条件稳定，但对非线性项仍然需要进行时间外插，因此，可将稳定外插方案SETTLS用于非线性项的计算。同时为了构建二阶精度的SISL时间积分方案，减小衰减作用，拟对线性项进行二阶精度非中央权重计算。

二时间层SISL方案的一般离散形式为

对非线性项采用SETTLS计算方案

Simmons等(1997)根据稳定性分析提出了一种线性项的二阶精度非中央权重计算方案，将该方案应用于线性项，则有

式中， 系数ξ决定衰减的大小，通常取ξ为0.05—0.1。可以证明该方案尽管对线性项取非中央权重，但仍具有二阶精度(Côté et al，1995)，且对低波数波动具有较小的衰减。

在采用半隐式时间积分方案时，首先须对模式方程使用参考大气廓线进行线性化。一般参考大气可有两种选择，一种为等温大气(T_r=常数)，另一种是等熵大气(θr=常数)，GRAPES模式取等温大气作为参考大气。

关于等温参考大气的选择，Simmons等(1997)指出，等温参考大气的温度要高于实际大气平均温度才能保证半隐式方案的稳定(中纬度地区的实际平均大气温度一般低于300 K)。Bénard(2003)利用简单的一维浅水波方程分析了二时间层SISL方案中参考大气与模式稳定性的关系，同样得到等温参考大气温度要高于实际平均大气温度才能保持稳定的结论。Temperton等(2001)将线性项的二阶精度非中央权重方案应用于欧洲中期天气预报中心(ECMWF)的二时间层全球预报模式中，其中，等温参考大气取T_r=300 K，ξ取为0.05。本研究参考ECMWF模式，在方案中采用Simmons等(1997)和Bénard(2003)的结论，取高于实际大气平均温度的等温参考大气温度(300 K)，以保证模式计算方案的稳定性。

Simmons等(1997)指出，线性项的一阶非中央权重系数与参考大气温度对稳定性具有等同的影响，增加非中央权重系数(即采用一阶非中央权重方案式(24))就不需要取较高的参考大气温度，这可能是GRAPES原SISL时间积分方案可以取较低的等温参考大气(T_r= 257.93 K)的原因。但是较大的非中央权重系数会对低波数波动产生较强的衰减作用。

由于使用SETTLS方案和新的SISL离散形式后上游点的计算方案和时间积分形式均发生了变化，因此，要对标量场和矢量的时间离散形式进行修改。将t时间层表示为n，t+Δt时间层表示为n+1，t+Δt/2时间层为n+1/2，标量场的SISL时间离散形式，对于热力学方程，可写为

式中，同样对于扰动无量纲气压Π′，则有式中，

计算上游点采用SETTLS方案后，不需要首先计算上游点与到达点的中间点，而是直接计算出上游点，因此需对矢量场的时间离散形式做一定的修改。设上游点D和到达点A其坐标分别为(λ_D，φ_D，z_D)和(λ_A，φ_A，z_A)，所对应的局地单位矢量分别为(i _*，j _*，k _*)和(i，j，k)，球面上D和A点的单位矢量之间的转换关系为

各个分量可写为 根据球面上单位矢量的转换关系则矢量场的时间离散形式变为式中， 式中， 式中，

值得注意的是，新SISL的时间离散形式仅改变了赫姆霍兹方程相应的系数，赫姆霍兹方程在形式上并没有发生变化。 4 理想试验 4.1 平衡流试验

对于SISL模式来说，拉格朗日轨迹的计算是整个模式动力框架的基础，关系到平流过程计算的合理性和正确性。由于本研究采用了新的上游点SETTLS计算方案，所以，需进行平衡流试验来验证上游点计算的正确性。参考薛纪善等(2008)设计的平衡流试验，初始水平风分量场u与气压场满足地转平衡关系，并呈平行的东西向的直线，详细的试验设计不再赘述。特别需要指出的是，试验初始场地面气温取为T₀=288 K，等温参考大气T_r=300 K，ξ取为0.05。

从新方案积分的结果(图 1)可以看到，模式积分15 d后水平风u(图 1b)仍然保持初始时刻理想场(图 1a)所具有的东西向平行结构，初始场的平行结构不随模式积分时间而改变。对于气压变量场Π同样保持初始场的平行结构而不随时间变化(图略)。平衡流试验的结果表明新的时间积分方案和拉格朗日轨迹计算是稳定和正确的。

从平衡流试验全区域内总动能((u²+v²+w²)/2)的时间序列(图 2)可以看出，原SISL时间积分方案的总动能随积分时间而减少，新SISL时间积分方案总能量虽也有很小的增加，但在16 d的积分时间内，总能量大体保持稳定。相比原SISL方案，新SISL时间积分方案总动能没有明显的衰减。可以看出新SISL积分方案有效地抑制了总能量的衰减。

4.2 非静力重力波试验

为了评估新SISL时间积分方案对重力波衰减作用的影响，设计了非静力重力波试验，采用孤立圆形山作为理想地形，地形廓线取为

式中，山体最大高度h_m=50 m，山体宽度半径为a=1000 m。模式水平分辨率Δx= 400 m，顶层高度为z_T=25.6 km，垂直分辨率Δ= 400 m，垂直取64层。假定初始时刻大气满足静力平衡关系，过山气流速度u为20 m/s，地面气温取为290 K，浮力频率N=0.01 s^-1，原SISL时间积分方案的参考大气温度T_r= 257.93 K，新SISL时间积分方案(式(26)、(27))参考大气取为等温大气T_r= 300 K。在模式顶部4 km对垂直速度进行隐式瑞利衰减(Klemp et al，2008)，消除重力波在刚性上边界反射对结果的影响。这时无量纲量Na/u=0.5＜1，根据Smith(1980)线性理论分析，气流经过山体时地形激发的重力波为非静力重力波。非静力重力波的结构特征是垂直倾斜向上传播，而静力重力波结构特征是垂直向上传播。

模式积分2 h后基本平衡，图 3给出了其重力波结构(垂直速度)分布。相比全球非静力模式如Model for Prediction Across Scales(MPAS)、Non-hydrostatic Icosahedral Model(NIM)(见Dynamical Core Model Intercomparison Project，DCMIP，https://earthsystemcog.org/projects/dcmip-2012/)模拟的非静力重力波结果，沿赤道一周大概有10个模左右闭合重力波数。本研究的GRAPES是有限区域模式，未使用周期边界条件，采用新SISL时间积分方案(式(26)、(27))后，闭合重力波数同样有10个左右(图 3b)，且在下游重力波呈明显的垂直倾斜结构，表现出了更多非静力重力波分量，这些特征均与全球非静力模式(如NIM，MPAS)的结果相近。而原GRAPES模式采用的SISL时间积分方案(式(24))，由于非中央权重系数ε=0.7，对重力波有较强的衰减作用，闭合重力波数明显减少，垂直倾斜度也明显减小(图 3a)，波动在下游被逐渐衰减，非静力重力波的分量并不明显，且总体上波动振幅明显弱于新SISL时间积分方案。

同样的结果亦可从垂直速度的水平分布上看出(图 4)，新SISL时间积分方案的重力波垂直速度扰动一直传播至格点340(图 4b)，而原SISL方案垂直速度的扰动仅传播至格点250附近(图 4a)。

5 总结与讨论

将两时间层稳定外插方案(SETTLS)引入GRAPES区域模式，并将其用于计算上游点和非线性项的时间外插计算，计算精度均为二阶。使用SETTLS方案计算上游点时，可直接计算出上游点的位置，而无需先计算出中间点。同时对于线性项采用了二阶精度的非中央权重计算方案，通过平衡流和非静力重力波等理想试验对新SISL时间积分方案进行评估，并与GRAPES原 SISL时间积分方案进行比较，得到如下结论：

(1)对线性项采用二阶精度的计算，要求等温参考大气温度必须大于实际大气温度，才能保证半隐式时间积分方案的稳定性。

(2)线性项的一阶非中央权重平均虽然对参考大气的温度没有要求(不要求取高于实际平均大气温度的参考大气)，但较大的非中央权重系数对低波数波动具有较强的衰减作用。对线性项采用二阶精度的非中央权重计算方案，同时取较高的参考大气温度(一般取为300 K)，使得线性项的计算具有二阶精度，对重力波的衰减明显减弱。

对于式(4)中ξ的选取，一系列的平衡流敏感性试验结果表明，在积分25 d以后选择较大的ξ能够增加模式的稳定性，但在短期天气预报中(一周以内)的选取对结果并无太大的影响，因此，对于一周以内的天气预报而言，倾向于取较小的ξ值，如0.05或0.1。

以上的结果均基于理想试验，在真实大气数值模式中尤其在应用于GRAPES区域模式后的预报效果需要做进一步的检验。