1 引言

自从Tapp and White(1976)提出完全弹性非静力模式(NH)以来，经过约30多年的努力，现在世界上的气象中心已建立的业务或准业务模式已不下10个，其积分范围大多是有限区域，少数既是全球又含有有限区域(廖洞贤，2004)。Steppler et al.(2003)、 Steppler(1995)和Dudhia(1995)认为: 垂直坐标的选择、 无摆动且计算精确的平流格式的构造和椭圆型方程的数值解的收敛速度应给予特别关注，其中第一个问题在采用高分辨率时因陡峭地形使坐标面倾斜会引起计算不稳定。不过，所谈内容仅限于计算方法，对于动力方面的问题还未涉及。

关于动力方面的问题，早期Tapp and White(1976)以及张可苏和周晓平(1980)在静力近似和完全弹性非静力情况下，比较了高频、 小振幅的重力惯性波的性质，结论是: 在静力情况下会引起虚假的高频重力波，应采用非静力方法。2002年，Yeh et al.(2002)还发现: 由于静力近似略去了dw/dt，非静力导致的内波崩溃和翻转作用当水平运动尺度在10 km以下时，不可忽视。但这些研究还限于较窄的范围，还没有进行一般性研究和能量分析。

近年来，一些气象中心企图用非静力模式制作全球中长期预报和气候预测，如UKNC(Cullen et al 1997)、 GEM(Cote et al., 1998)和GRAPES(薛纪善等，2008)(缩写词全称见附录，下面其他缩写词全称亦然)。它们的预报效果表明: 与采用同样分辨率的静力原始方程模式(全球谱模式或差分模式)相比，目前还看不出有明显的优势(Yeh et al，2002)。另外，用GRAPES模式计算的垂直速度有时会异常大，引起虚假的格点暴雨。WRF模式也反映有类似问题(Skamarock et al., 2005)，WRF曾用约束Courant数以达到约束垂直速度的nudging方法。

不论对于静力模式或非静力模式，能量守恒和能量转换都应是模式设计中需要考虑的。对此，很多人都进行过研究，特别是对于前一种模式。由于这类模式发展的时间较长，结果比较成熟，可以作为后一类模式设计的参考。早在20世纪60年代，Arakawa(1966)就提出过正压过滤模式长期预报积分的方案，以后，Arakawa and Lamb(1977)、 Arakawa and Suartz(1983)又研究UCLA大气环流模式的保持能量守恒的数值方案和sigma坐标系中原始方程模式的垂直差分方案，其结果多次为数值预报模式设计者所引用。但是，Arakawa等的方案只能在时间连续，空间离散情况下保持能量守恒(即瞬时能量守恒)，这是其不足之处。对于这个问题，曾庆存和张学洪(1987)以及左瑞亭等(2004)曾提出时空差分的能量守恒的隐式方案。1995年，季仲贞等(1995)提出时空差分的能量守恒的显式方案。至于完全弹性非静力模式，由于问题比较复杂，研究能量守恒的较少。研究方案有针对有限区模式垂直差分的Taylor(1984)方案和针对有限区模式的Gallus and Rancic(1996)方案，以及致力于全球气候模式的Satoh(2002)方案。Satoh以连续方程、 动量方程和全能量方程构成模式的动力方案，计算结果令人鼓舞，其质量守恒和总能量守恒的计算都在舍入误差精度以内。

根据上面的分析和提出的问题，下面我们将作进一步的分析并对某些问题提出解决的方法。全文共6节，第2到第6节依次是: 对垂直运动方程进行尺度分析的结果和垂直扰动运动方程各项量级接近的可能、 垂直扰动气压梯度力的截断误差、 多级扰动法、 多级扰动的性质和结论。

2 对垂直运动方程进行尺度分析的结果和方程各项量级接近的可能

根据廖洞贤(2008)研究中的尺度分析结果，我们可以把在绝热无耗散情况下大气垂直运动方程各项的量级表示出来(图 1)。图 1中纵坐标是lg10a，10a是量级，横坐标是方程各项。从图 1可以看出如下几个特点:

(1)不论对于水平尺度L=10⁶ m、 L=10⁵ m或L=10⁴ m的运动，最大项都是垂直气压梯度力项(VPGF)和重力项(G)，而最小项却各有不同，对于L=10⁶m和L=10⁵m的运动，是随L的减小而增大的垂直平流项(VAD); 对于L=10⁴m的运动，是曲率项(CUR)。由于VPGF、 G和CUR不随L而变，VPGF 与CUR之差也不随L而变，故对于这三种水平尺度的运动，最大项和最小项之间的量级差(以下简称最大量级差)总是随L的减小而减小的。根据图 1所示，对于L=10⁶m、 L=10⁵ m 或L=10⁴ m的运动，最大量级差依次是9个、 7个和6个量级。

(2)最大项和次大项之间差4个量级。根据廖洞贤(2008)的研究，VPGF和最小项之间的量级差最大，应取的垂直分辨率最高，而对于一定的L，带有垂直差商项的量级愈大，其截断误差的量级也愈大。因此，通过分离技术使VPGF相应扰动项的量级减小，会使所需垂直分辨率降低，而在同样分辨率下，用分离技术时的截断误差比未用时的截断误差显著减小。在已设计的完全弹性非静力模式中，Tapp and White(1976)、 Dudhia(1993)、 Bubnova et al(1995)、 Xue et al.(2000)和薛纪善等(2008)的模式，毫无例外地都采用只随高度而变的气候值作基本态变量，并用来对随时空而变的热力学变量进行分离，用对后者的扰动作预报。看来这样做有一定好处。为了了解这种好处，廖洞贤和孙岚(2003)曾用19种地面温度(T_s)和温度的垂直递减率(Γ)的组合构成的温度廓线和用气候值且只是z的函数的基本态廓线所作的计算，扰动项的量级比VPGF和G仅降低约1～2个量级，还不足以使方程达到各项量级接近的程度。但是，也应看到如果我们用类似方法对扰动进行分离，则有可能使更高一级扰动项的量级继续降低，经过几次使方程达到各项量级接近，使与最大量级差相应的垂直分辨率达到可用的程度。当然，这些想法还须分析论证。

3 垂直扰动气压梯度力项的截断误差3.1 离散的垂直扰动运动方程和其各项的量级

设δF/δz=əF/əz+E(δF/δz)，其中E(δF/δz)是δF/δz的截断误差。以L=10⁶ m的运动为例，我们可以把离散的大气垂直扰动运动方程写成:

这里，方程中各项下面的数表示各项的量级 [参阅廖洞贤(2008)的方程(18)及廖洞贤和孙岚(2002)]。 F⁽¹⁾是热力学变量F对基本态变量 的1级扰动，即F= +F⁽¹⁾。因为截断误差的量级比待离散量的量级为小，我们不妨取水平截断误差项E_H为方程的水平待离散项的1/10，这样E_H为10^-8。E_V是方 程中垂直截断误差的总和，它可以用垂直扰动气压梯度力(VPPGF)项的截断误差E(VPPGF)近似地表示。

3.2 E(VPPGF)对dw/dt或əw/ət的影响

在方程(1)中，除被离散项外，如果不计观测、 分析或舍入误差的影响，其他项都是可以精确计算的。这样，在第1节中提到的 w异常增大应该是dw/dt或əw/ət出现异常造成的，而dw/dt或əw/ət出现异常只可能是方程中出现很大 E(VPP-GF)所致。如果E(VPPGF)是VPPGF的1/100，则E(VPPGF)的量级是10^-2。于是，E(VPPGF)的量级就大于fu、(u²+v²)/a和平流项的量级。

3.3 对E(VPPGF)的估计

为了得到具体的结果，我们还须对E(VPP-GF)进行估计。 取T和静止大气基本态变量各是

和

其中，下标s和t表示地面和模式大气顶; 变量上带^的量表示基本态变量，不带^的量表示实际变量; 下标T1、 T2各表示实际变量和基本态变量在对流层顶的值; z₁表示T廓线中的等温层顶。取z_s=0 m，z₁=3×10³ m，z_T1=1×10⁴ m，z_T2=1.2×10⁴m，T_s=278.15 K，Γ₁=5.5×10^-3Km^-1， _s=288.15 K，Γ₂=6.5×10^－3 K/m。于是，从静力方程我们有

和

其中，，是气压在z₁处的值; 同样，p_T1和分别表示在z_T1处的气压和在z_T2处的基本态气压。这样，我们可以求出p⁽¹⁾，ρ⁽¹⁾。至于E(VPPGF)，在取中心差的情形，利用Taylor展开可以得到

其中，k是垂直网格的序号。于是，从式(6)，可以求得在Δz=1000 m，200 m和20 m时的|E(VPPGF)| ，它们各可以表示其E(VPPGF)的量级，如表 1所示。

比较Δz=1000 m时的|E(VPPGF)|和方程(1)中各项的量级，可以看出: 除VPPGF和浮力项(BUOY)外，它比所有的项都大。这样，əw/ət 受到E(VPPGF)的影响就掩盖了fu，(u²+v²)/a和平流项等对əw/ət 的贡献。而且，由于VPPGF总是方程中的最大项，这种现象还会持续存在; 虽然(u²+v²)/a总为正，fu在北半球西风带中基本为正，它们的持续作用仍不容易表现出来。这是不合理的。

3.4 垂直分辨率的问题

从前面的分析讨论来看，E(VPPGF)所以不合理，是选取垂直分辨率不当造成的。其实并不尽然。从表 1可以看到: 提高垂直分辨率可以使E(VPPGF)减小，但是，取到Δz=20 m仍有5.0384×10^-6 m/s²的误差; 对于L=10⁶ m的运动，对əw/ət的影响仍大于水平和垂直平流项对əw/ət 的贡献。 看来要彻底避免这种现象的出现，只有再提高分辨率才行。不过，根据朱禾等(2007)，要避免这种现象，Δz须取为

其中，H是大气垂直厚度的特征尺度，取为10⁴ m，Az和B₁各是VPPGF和最小项的量级(10^a)的次数，(Az－B₁)M=max (Az－B₁)_，是方程的序号，r是调节精度的参数，取为1。在非静力模式中，(Az－B₁)M就是方程(1)中的Az－B₁即8。根据式(7)的计算，对于L=10⁶ m和L=10⁵ m的运动，Δz须取到 0.77 m和7.75 m。这显然在当前客观条件下是很难实现的。而且，即使实现，取这样小的Δz，舍入误差反而会增大，使计算效果适得其反(冯康等，1978)。不过还可以采用第4节中的多级扰动法。

4 多级扰动法 4.1 多级扰动的构造

根据方程(1)，和是其最大项。在这个条件下，我们可以得到如下结果:

结果1 和量级相等符号相反。

因为，如二者量级不相等，则其中必有一项的量级大于另一项，这样方程中最大项只有一项，这就违反了尺度分析中一个方程的最大项至少有2项的结论(叶笃正等，1988; 刘式适和刘适达，1991); 另外，如它们的符号相同，则方程不能达到平衡，故结果1成立。

根据结果1，我们可以构造1级扰动平衡方程 

其中，p⁽¹⁾，ρ⁽¹⁾各是p⁽¹⁾和ρ⁽¹⁾中满足上方程的部分，p⁽¹⁾=p⁽¹⁾+p⁽²⁾，ρ⁽¹⁾=ρ⁽¹⁾+ρ⁽²⁾。但是 ⁽¹⁾和⁽¹⁾还必须各和，一起满足形如

的状态方程才能和θ₁保持一致。这里p₁=+p⁽¹⁾，ρ₁=+ρ⁽¹⁾，θ₁=+θ⁽¹⁾，pr是参考气压，取为常数。

方程(8)推广，设

其中，n表示扰动级次，n是大于或等于1的正整数; F⁽ⁿ⁾表示F⁽ⁿ⁾中满足n级扰动平衡方程

又和p_n－1、 ρ_n－1、 θ_n－1一起满足状态方程

的部分，其中

其中n是大于等于0的整数; 当n=0时p₀=，ρ₀=，θ₀=。

而F⁽ⁿ⁺¹⁾表示F⁽ⁿ⁾对F⁽ⁿ⁾的扰动，如取

其中c₁是一近于1的正实数，0.9≤c₁≤0.99，则

而从方程(11)，如ρ_n、 p_n和、0⁽¹⁾，直到0^(n－1)已知，还可以求出0⁽ⁿ⁾，从而，从ρ⁽ⁿ⁾、 θ⁽ⁿ⁾、 p⁽ⁿ⁾通过(9)求出ρ⁽ⁿ⁺¹⁾、 θ⁽ⁿ⁺¹⁾、 p⁽ⁿ⁺¹⁾。不过，在工作中还可以利用一些近似式来达到目的。在下节中我们将进行推导。

4.2 一些近似关系

对状态方程

取对数，有

. 而对方程(11)，置n=1，取对数，并利用上式，则有

故近似地有

， 或

其中 。类似地，用p、 ρ、 θ代替方程(11)中的p_n、 ρ_n、 θ_n利用p、、 p⁽¹⁾等的关系，我们还有

从方程(11)取对数，考虑其线性部分，还有

其中，。如p=p_n+p⁽ⁿ⁺¹⁾，ρ=ρ_n+ρ⁽ⁿ⁺¹⁾，θ=θ_n+θ⁽ⁿ⁺¹⁾，还有

故

其中(c²_s)_n=γR_dT_n。

4.3 扰动项量级随级次的变化

用ρ乘方程(1)，得到

利用式(9)和方程(10)，置n=1，则方程(21)成为

易见，əp⁽²⁾/əz和ρ⁽²⁾g仍是方程(22)中仅有的最大项。于是根据推导结果1同样的理由，我们可以推得: əp⁽²⁾/əz和ρ⁽²⁾g量级相等，符号相反的结论。于是，我们可以构造2级扰动平衡方程 

而

故

这里<_M表示量级小于; ～表示量级相等。这样一步一步地推下去，当n>2时，如扰动方程是

则只要ρ⁽ⁿ⁾g>_Mρfu，ρ⁽ⁿ⁾g和əp⁽ⁿ⁾/əz都是方程中仅有的最大项。于是结果1可以推广为:

结果2 əp⁽ⁿ⁾/əz和ρ⁽ⁿ⁾g量级相等，符号相反。

于是，利用推导不等式(24)和(25)同样的方法，可以得到

且

式(27)、(28)说明: 随着n的增加，扰动项量级随着级次的增加而减小。但是，扰动方程中当原方程(1)中的次大项fu成为n级扰动方程的最大项时，方程中的最大项可能不只2项，而是3项，这时结果2不一定成立。为了方便，不妨设结果2不一定成立的n是n₁。

为了说明多级扰动法的作用，我们利用扰动项量级随n减小的性质来估计E(VPPGF)。用n代替式(6)中的1，用3.3节中的理想场，按式(9)、(10)、(12)、(13)和附录中解n级扰动平衡方程的方法，仍取表 1中的分辨率，可以求出|E(VPPGF)|，结果如图 2所示。从图 2可见，|E(VPPGF)| 的量级随Δz的减小而减小，也随n的增加而减小，都可以使E(VPPGF)对倾向的影响减小。容易看出: 不论取c₁=0.91或0.95，当n=3，Δz=200 m时，|E(VPPGF)| 都小于10^-5，|E(VPPGF)| 对倾向的影响已不致歪曲或掩盖fu和(u²+v²)/a对倾向的贡献了。

4.4 扰动项量级减小和垂直分辨率的关系

如果n增加1，扰动项减小一个量级，则扰动项ρ⁽¹⁾g/ρ减小到和u的量级相等时，fn=n₁，这时VPPGF减小了3个量级。按式(7)用2阶中心差应取的Δz如表 2所示。如果认为Δz大于100 m为可用，则从表 2可见，除L=10⁶ m的运动外，对于L=10⁵ m和L=10⁴ m的运动，Δz可用。

至于L=10⁶ m的运动，如我们只要求E(VPPGF)的量级小于曲率项，则VPPGF和BOUY只减小2个量级，(Az－B₁)M为3，按式(7)，Δz= 245 m，也是可用的。

5 多级扰动的性质 5.1 2 个关于 Δ(=p⁽¹⁾/z+ρ⁽¹⁾g)的结果

结果3 利用式(9)和方程(10)有

上式说明: 在应用多级扰动法时，尽管n可以变化，但与对流活动密切相关的Δ不随n而变，垂直扰动运动方程中其他项也不随之而变，故n的变化不改变方程的性质。

结果4 设有满足静力方程的变量p₁、 ρ₁和p₂、 ρ₂，则它们的差p′(=p₁-p₂)和ρ′(=ρ₁-ρ₂)必满足式 

易见，如p₂、 ρ₂为只随z而变的静止大气基本态变量，则有

5.2 静力扰动 5.2.1 静力扰动的性质

显然扰动平衡方程中的气压扰动 ⁽ⁿ⁾是静力扰动。以n级扰动平衡方程 

为例，易见，在的条件下(证明见5.3节)，可以得到

其中，J表示一个很高的气层。

可见，p⁽ⁿ⁾表示以| p⁽ⁿ⁾|为密度的静止大气单位截面，从高度z到∞的空气柱的正(负)重量扰动，即正(负)静力气压扰动。故 ⁽ⁿ⁾是静力扰动。易见，这个结果可以推广为: 方程 的解p⁽ⁿ⁾也是静力扰动。

5.2.2 多级扰动静力扣除法

多级扰动法实质上是多次利用扰动平衡方程的解，从p⁽ⁿ⁾中扣除其所含大小与p⁽ⁿ⁾相同的静力扰动部分，使p⁽ⁿ⁾量级减小，达到əp⁽ⁿ⁾/əz与其它项量级接近的一种方法。这种方法也可以称作 “多级扰动静力扣除法”。

5.2.3 修正的大气基本态变量

从大气总的来看，因为ρ=ρ_n+ρ⁽ⁿ⁺¹⁾，p=p_n+p⁽ⁿ⁺¹⁾，易证əp_n/əz=－ρ_ng，故ρ_n、 p_n各可代替和作为大气修正的基本态变量，ρ⁽ⁿ⁺¹⁾和p⁽ⁿ⁺¹⁾各代替ρ⁽¹⁾和p⁽¹⁾作为ρ、 p和ρ_n和p_n之间的扰动。这种基本态变量虽然不一定是z的函数而是时间和空间的函数，但因ρ⁽ⁿ⁺¹⁾=(1-c₁)ⁿρ⁽¹⁾，p⁽ⁿ⁺¹⁾=(1-c₁)ⁿp⁽¹⁾使ρ⁽ⁿ⁺¹⁾ρ⁽¹⁾，p⁽ⁿ⁺¹⁾p⁽¹⁾，当n≥1时，这种代替使ρ_n、 p_n大气垂直廓线和ρ、 p大气垂直廓线或实际大气垂直廓线达到很接近的程度; 而且，每次修正都是通过静力平衡达到的，可以认为: 模式大气的绝大部分已不含有垂直声波。不过，由于(u²+v²)/a，u和dw/dt的存在，它们还会引起非静力扰动，从而，使这种作用降低。不过，非静力扰动一般远比静力扰动为小，其所起的作用一般不太明显。只是，这时仍存在水平声波，考虑计算稳定性问题，应当把它作为主要的抑制对象。下面我们用一个例子来说明。

5.2.4 实例计算

为了说明修正的大气基本态变量垂直廓线和实际大气变量垂直廓线当n≥1时相接近的程度，我们把NCEP的2003年7月31日1800 UTC在(30°N，120°E)从1000 hPa～10 hPa的等压面温度资料，通过线性插值转换为z坐标的以200 m为Δz的格点值T，p按Lorentz网格依次放置在整数层和半数层上。还采用3.3节用理想场计算的z坐标格点的和，用前面的方法算出ρ_n和p_n，从而算出ρ－ρ_n和p－p_n。但这时除ρ－ρ_n是迅速随n减小外，p－p_n从n=0到1减小得很快，而从n=1到3则随n有所增加，增加值占从n=0到1减小值的15.8%(图略)。其所以如此，显然和计算中求z坐标p、 ρ值时采用的线性插值有关，因为这种插值不一定满足静力平衡条件，即使用原NCEP温度资料得出的ρ、 p格点值满足该条件。由于这样，我们不能用公式p⁽ⁿ⁾=(1-c₁)^n－1p⁽¹⁾来估计p－p_n，而只能用4.1节中p_n的表示式计算，这样，在p－p_n中不仅含有插值误差，还含有因计算 ⁽¹⁾等引入的数值积分的误差，这比ρ－ρ_n中所含的误差还多。为了这个缘故，我们改变了计算方案: 只对ρ插值，而p则用静力方程求出，并令其他方法不变，和也不变。这样得到的ρ－ρ_n和p－p_n的结果各表示在图 3和图 4中。从这两幅图中可以看出: ρ－ρ_n和p－p_n都随n减小，只是当n从0到1时减小得很快，而1到3时，减小的幅度较小罢了。这说明: 修正的大气基本态变量垂直廓线和实际大气变量垂直廓线是随n的增加而接近的。注意: 为了使单位协调，ρ采用hPa · s² · m^-2，p采用hPa(Stull，1988)。

5.2.5 非静力扰动

凡形式上不满足形如方程(32)的扰动都是非静力扰动。如

为了方便，我们不妨把式(32)改写成

并令p⁽ⁿ⁾=p⁽ⁿ⁾H+p⁽ⁿ⁾N，其中下标H，N各表示静力变量和非静力变量; S⁽ⁿ⁾是z的函数; 且 

把方程(36)代入方程(35)，则 

设，则

这是S⁽ⁿ⁾引起的非静力扰动，它是对流发展的重要因素。

5.3 扰动在∞高处的极限

设p，和ρ，均各满足5.1节中结果4和Δ=0的条件，则可以推出如下结果。

结果5 。其中 ; 。证明如下。先证n=1的情形。这时

由于在大气中1/T(ξ)有界，且总大于零，故

同理

且

故

因 

故式

且

故式(43)在区间 ［0，∞)成立。至于p⁽²⁾、p⁽²⁾，因p⁽²⁾=p⁽¹⁾－p⁽¹⁾，故

故用和前面同样的方法可以证式p⁽²⁾=c₁p⁽²⁾在区间 ［0，∞)成立。如此继续可证

和

结果6 。证明如下。对方程(42)作从zt到∞处的积分，利用式(44)，有

又因从n级扰动平衡方程

和式(48)，有

注意计算p⁽¹⁾t和⁽ⁿ⁾t时须化为数值积分的形式，且积分上限还须用有限高度代替。

6 结论

从前面的讨论，可以得出以下初步结论:

(1)用多级扰动法可以使大气垂直扰动运动方程中量级最大的垂直扰动气压梯度力项和浮力项的量级随级次减小，使方程中最大的截断误差对倾向的影响不致歪曲或掩盖较小的垂直科氏力项和曲率项对倾向的贡献，从而使倾向计算合理，可以避免w出现异常增大的现象; 而具有持续作用的垂直科氏力项和曲率项的作用也可以得到充分发挥，有利于制作较长时期的天气预报。

(2)扰动平衡方程的解表示静力扰动。多级扰动法实质上是多次利用不同级次的静力扰动扣除原扰动中和p⁽ⁿ⁾大小相同的静力部分，达到扰动方程中各项量级接近的一种方法。而且，从大气总的情况来看，经过这法不多几次修正，大气基本态ρ_n，p_n垂直廓线的绝大部分和实际大气ρ，p垂直廓线已很接近，可以认为，通过静力平衡，模式大气基本上已不含有垂直声波; 抑制快波使计算稳定的主要对象应是水平声波。

(3)可以证明: 在应用多级扰动法时，尽管级次n可以变化，但与对流活动密切有关的Δ(=əp⁽ⁿ⁾/əz+ρ⁽ⁿ⁾g)却不随n而变，因各级扰动垂直运动方程中其他项均不随n而变，故n的变化不改变原扰动方程的性质。

附录 1 n 级扰动平衡方程的数值解

n级扰动平衡方程可以写成 

对方程(A1)作从z=0到z_t的积分，采用Lorentz网格，可以得到

再作从z=z_{k－ 1/2} 到z_t的积分，可以得到

于是，可取

从式(13)和式(51)有

和

2 高级扰动的构造

(1)Δ=0的情形。根据式(47)，从式(A5)和(A6)可以得到

和

至于ρ⁽ⁿ⁾_k、 ρ⁽ⁿ⁾s和θ⁽ⁿ⁾_k、 θ⁽ⁿ⁾s等，则可以从式(13)和4.2节中有关公式推出。下同。

(2)Δ≠0的情形。这时从式(9)有

3 缩写词的全称

缩写词的全称如表 A1所示。