早在20世纪90年代末，数值预报专家就提出了物理过程和动力框架之间耦合的重要性，指出物理过程和动力框架之间的耦合对预报的精确性和准确性有至关重要的影响^[1-4]。这些研究细致分析了物理过程和动力框架之间耦合在时间上的协调性，探讨半隐式半拉格朗日大步长积分方案下物理过程和动力框架如何耦合才能得到准确解，提高耦合精度，从而改善模式的综合预报效果。

不仅是在时间上，在空间上数值预报模式也存在物理过程和动力框架耦合的协调性问题。目前，绝大多数数值预报模式的物理过程虽然在水平方向上单点描述，但垂直方向格点离散化方案对物理过程参数化表述方法和耦合有很大影响。数值模式垂直方向通常有两种离散方案：一种是Lorenz跳点，与物理过程有关的量，如温度、湿度、水物质、水平风等均分布在半层^[5]；另一种是Charney-Phillips跳点(简称C-P跳点)，水平风分布在半层，温度、湿度和水物质分布在整层^[6]。研究表明:在强涡流和强层流系统、地转适应过程等问题的处理中，C-P跳点优于Lorenz跳点设计^[7]。目前国际上大多数气象业务和科研中心的数值天气预报模式采用Lorenz跳点，因此, 多数物理过程根据Lorenz跳点分布设计。但也有少数的气象业务中心数值天气预报模式采用C-P跳点，如英国气象局。

GRAPES模式^[8-15]的动力框架考虑到C-P跳点的优势，因此，在垂直格点设计上采用C-P跳点设置，所有物理过程^[16]均在Lorenz跳点上计算。鉴于这样的设计，为了解决动力框架和物理过程在垂直格点分布上的不一致，GRAPES模式在进入物理过程之前，将位温和水物质从整层插值到半层；在物理过程计算完毕之后，将物理过程反馈给动力框架的位温和水物质同样再从半层插值回整层。尽管在插值过程中采用高精度的插值方案并考虑插值方案的单调性，在一定程度上减少了插值误差，提高了模式的预报效果^[17]，但物理过程调用前后的插值无法完全避免插值误差，同时由于所用插值方案及边界层处理无法保证边界层参数化方程中二阶导数的光滑性，造成物理量分布的空间不连续性。因此, 对于边界层过程，需要考虑设计在C-P跳点的参数化方案解决其与动力框架之间的耦合，以避免插值造成的误差和物理-动力耦合的不协调。

由于C-P跳点的温度和水平风是跳层分布，使湍流通量和地表通量物理过程与动力框架之间的耦合相对复杂^[14]，边界层参数化方案中直接考虑湍流通量在C-P跳点的表述也会遇到诸多问题。Holdaway等^[18-19]提出了一系列稳定边界层中局地闭合在C-P跳点下的可能参数化方法，详细比较了Lorenz跳点和C-P跳点下边界层过程和动力框架耦合的空间协调性，讨论了边界层在哪种格点上进行参数化更占优势。对于GRAPES模式中的C-P跳点，需要研究边界层过程在C-P跳点上的描述方法，避免插值误差的产生，解决边界层过程和模式框架耦合的协调性。

本研究将基于GRAPES全球模式现有的边界层方案，分别讨论稳定边界层和局地不稳定边界层的C-P方案设计以及边界条件问题，重新设计C-P跳点上边界层参数化方案，以解决边界层过程和动力框架垂直格点不一致的问题，使边界层和动力框架之间的耦合相协调。通过本研究可初步了解物理过程调用前后两次插值对模式性能的影响。

大气边界层(PBL)过程是大气运动物理过程中的一个重要组成部分。大气运动方程组由于Reynold分解和平均产生了u′w′, w′θ′等非线性未知量，其中, u′为水平风扰动量，w′为垂直气流扰动量，θ′为位温扰动量，u′w′, w′θ′代表湍流运动对动量、热量的垂直输送。因为湍流的尺度较小，数值模式的网格无法显式描述这些新出现的未知变量，所以大气边界层参数化要使用平均量描述这些次网格的通量，使方程闭合，从而得以求解。这种用平均量描述这些次网格边界层通量的过程叫做边界层参数化。

目前GRAPES全球模式边界层方案采用的是新的MRF边界层方案，该方案考虑一阶的K闭合条件，使用K廓线闭合描述不稳定边界层中的非局地扩散过程^[20]，同时方案考虑了非局地扩散造成的反梯度输送的影响^[21]；在稳定大气中，方案采用局地K闭合描述局地不稳定扩散过程。美国NCEP中期预报模式中边界层过程的描述就是基于该方案^[22]，即新的MRF方案。在新的MRF方案中考虑了边界层顶云的夹卷效应对湍流输送的作用，使有云情况下的湍流输送描述更为合理^[23-24]。

边界层参数化方程为

(1)

式(1)中, ψ表示动量、温度、湿度以及云水和云冰, t表示时间, w′ψ′表示湍流运动对动量、热量等的垂直输送, z表示垂直高度, K表示热动量扩散系数, 反梯度输送项γ表示大尺度涡对于通量的贡献。求解式(1)的关键是求取相应的交换系数K，然后对式(1)进行垂直格点离散化，隐式求解即可得到边界层物理过程新一时刻的ψ，从而得到边界层物理过程的倾向。

表 1给出了GRAPES全球模式边界层方案中各变量在C-P和Lorenz垂直格点上的分布。其中θ为位温，q为比湿，u和v为水平风，K_H为热量交换系数，K_M为动量交换系数，(w′θ′)_s为地表热通量，n为垂直层次，N为顶层。当在C-P跳点上设计边界层方案时，水平风的离散化方程不变。温湿的离散化方程在整层上求解，其交换系数K_H必需在半层上求解。

表 1

表 1 边界层方案变量在C-P和Lorenz垂直格点上的分布 Table 1 C-P and Lorenz configurations of variables in PBL scheme 垂直层次 C-P Lorenz n=N(上边界) θ, q n=N-0.5 u, v, KH θ, q, u, v ┋ ┋ ┋ n=1.5 u, v, KH θ, q, u, v n=1 θ, q，KM KH, KM n=0.5 u, v, (w′θ′)sKH θ, q, u, v n=0(地表) θ, q, u* w′θ′s, u*

表 1 边界层方案变量在C-P和Lorenz垂直格点上的分布 Table 1 C-P and Lorenz configurations of variables in PBL scheme

如表 1所示，C-P跳点相较于Lorenz跳点的区别在于温度和水平风是交错格点分布的，因此，热量交换系数K_H和动量交换系数K_M也是交错格点分布。C-P跳点的地表热通量根据常值通量层理论在n=0.5半层上分布。

模式中近地层到边界层顶就是充分混合层，这一范围采用非局地K廓线闭合方案求解扩散系数。动量扩散系数K_M和热量扩散系数K_H分别为

(2)

(3)

其中，k=0.4为冯卡门常数，w_s为近地层速度尺度，z为垂直高度，H为边界层高度，P_r为Prandtl数，不稳定混合层中P_r=Φ_H/Φ_M，Φ_H和Φ_M分别为热量和动量的普适函数。

由式(2)和式(3)可以看出，不稳定边界层中的交换系数仅与近地层的不稳定度和边界层高度有关，因此，在C-P跳点上，不稳定边界层中的热量交换系数K_H只需将高度z由整层更换到半层即可得到半层的热量交换系数，而动量交换系数K_M无需改变。

由于C-P跳点的水平风和温度不在同一个垂直层次，所求的热量和动量交换系数也不在同一层次，因此，稳定条件下如何计算局地理查森数R_i，如何求取交换系数K_H成为问题的关键。

在Lorenz跳点上，通过半层的虚位温θ_v，就可以求出整层的浮力，从而得到整层上的局地理查森数R_i：

(4)

式(4)中，g为重力加速度，U为全风速。注意R_i与水平风和温度都有关。整层上的交换系数为

(5)

式(5)中，l为混合长，F_{H, M}是R_i的经验函数。

在C-P跳点上水平风和温度不在同一层，因此，要计算某一层的R_i，或将与风相关的量插值，如风剪切，得到整层的理查森数；或将与温度相关的量插值，如浮力，得到半层的理查森数；或将两者分别插值，得到整层和半层的理查森数，从而分别得到动量和热量的交换系数。

根据已有的经验^[18-19]，将温度相关量插值求取理查森数的效果最好。因此，首先用整层的温度计算出半层上的浮力：

(6)

再通过线性插值得到整层上的浮力B，求取整层上的热量交换系数和动量的交换系数K_M。再次通过线性插值，求得半层上的热量的交换系数K_H。

下边界条件是边界层方程求解的关键，直接影响由地表驱动的从下向上的湍流交换。Lorenz跳点中边界层方程在半层上求解，其下边界条件直接用地表通量给出。而在C-P跳点上，温湿方程在整层上求解。利用常值通量层的假设，假定从地面到模式最低层(半层)上的通量不变，温湿方程的下边界条件同理近似用地表通量给定。该假定在模式最低层足够低的情况下是可以满足的，一般认为模式最低层应小于边界层高度的十分之一。目前GRAPES第1层整层为22 m，模式最低层即为11 m，上述假定可以成立。

由于温湿的离散化方程是从地表向上一层开始建立，地表层上不求解方程，直接将地表上一层的温湿倾向值赋给地表层。这在近地表处垂直分辨率较高的情况下也是合理的。

将C-P边界层方案模块接到GRAPES_GFS中，进行了个例分析，并与原方案在Lorenz跳点调用的结果进行比较。

选取NCEP FNL分析资料作为初始场，起报时间为2013年5月1日12:00(世界时, 下同)。模式分辨率为0.25°×0.25°，垂直层次60层，时间步长取300 s。物理过程选项包括RRTM辐射方案、CoLM陆面模式、新的SAS积云对流参数化方案，重力波拖曳使用欧洲中期天气预报中心方案，微物理为GRAPES自主研发的复杂冰相方案。控制试验为所有物理过程均在Lorenz跳点上调用，边界层用NMRF方案，对比试验边界层用C-P跳点方案，且在C-P跳点上调用。

图 1显示了预报6 h的全球平均温度和湿度的边界层倾向。由图 1可见，边界层在Lorenz跳点上调用由于插值导致不稳定边界层内的虚假振荡非常明显，尤其是湿度倾向的振荡尤为强烈；而C-P边界层方案的温湿倾向非常光滑，说明边界层过程在C-P跳点上计算和调用有效地抑制了这一不合理的振荡。

图1

图1 预报6 h的全球平均边界层过程瞬时位温和湿度倾向廓线(a)位温倾向廓线，(b)湿度倾向廓线 Fig.1 Global averaged instantaneous potential temperature and moisture tendency profiles for planetary boundary layer (PBL) at 6 h forecast (a) potential temperature tendeny, (b) moisture tendency

图 2给出的是6 h边界层瞬时湿度倾向的纬向平均，可以看出, 控制试验中边界层方案有着非常明显的条纹，显然是不合理的，而C-P边界层方案很好地消除了这一条纹噪音，使边界层湿度倾向分布更为合理。图 2b中能清楚地看出热带地区湿度倾向呈双峰结构，这正是边界层过程中地表驱动和边界层顶积云驱动共同作用产生的湍流交换造成的。

图2

图2 预报6 h的边界层过程瞬时湿度倾向纬向平均(单位：10-4 g·kg-1·s-1) (a)控制试验，(b) C-P边界层方案 Fig.2 Zonal mean instantaneous moisture tendency for PBL at 6 h forecast (unit:10-4 g·kg-1·s-1) (a) CTRL experiment, (b) C-P PBL scheme

图 3显示的是预报6 h全球平均的温度场和湿度场，发现原来的控制试验不仅是温湿倾向在低层有强烈的虚假振荡，温湿廓线本身在低层也存在振荡的不连续分布，尤其是湿度场，可以看到, 低层的振荡非常明显，模式第3层的湿度甚至比第4层还要低，这显然是不合理的。而应用了C-P边界层方案之后，温度场和湿度场的振荡消失，廓线光滑很多，其结果更为合理。

图3

图3 预报6 h的全球平均温度和湿度廓线(a)位温廓线，(b)湿度廓线 Fig.3 Global averaged potential temperature and moisture profiles at 6 h forecast (a) potential temperature, (b) moisture

从很多单点分析可以更清楚地看出, 控制试验中温湿廓线存在的锯齿状振荡以及C-P边界层方案改进的效果。图 4给出了30°S，60°E的温湿廓线，可以看出, 该处边界层高度约为1000~1500 m，控制试验中该点的温湿廓线在不稳定边界层中有非常明显的锯齿状振荡。

图4

图4 预报6 h的30°S，60°E点的位温和湿度廓线(a)位温廓线，(b)湿度廓线 Fig.4 Potential temperature and moisture profiles at 6 h forecast at 30°S, 60°E (a) potential temperature, (b) moisture

大气运动具有高度的非线性特征，大气模式系统本身具有初值敏感性^[25]，模式预报中边界层温湿廓线的锯齿状振荡相当于在模式低层不断引入噪声扰动，这对于中期预报的结果具有非常不利的影响。而边界层过程在C-P跳点计算和调用会有效地抑制该振荡，垂直廓线要光滑很多，其结果更为合理。

进行了15 d的连续试验，参数设置同上，时间是2013年5月1日12:00-15日12:00，预报时效为8 d。通过统计检验发现C-P边界层方案对于全球、北半球及东亚的中低层形势场均有一定改进。由于C-P边界层方案直接影响的是低层的温湿场，首先分析其对低层温度的影响。图 5显示的是北半球850 hPa温度场的均方根误差。可以看出，C-P边界层方案对低层的温度场稍有改善，前3 d的改善达到0.05显著性水平。

图5

图5 北半球850 hPa温度场的均方根误差检验(a)均方根误差，(b)显著性水平 Fig.5 The root mean square error of 850 hPa temperature in the North Hemisphere during 1-15 May 2013 (a) root mean square error, (b) significance level

图 6显示了北半球地区500 hPa温度场的距平相关系数的比较，由图 6可见，C-P边界层方案的应用在北半球地区对500 hPa温度场也有一定改进。且前4 d的改善达到0.05显著性水平。

图6

图6 东亚地区500 hPa温度场距平相关系数对比检验(a)距平相关系数，(b)显著性水平 Fig.6 The anomaly correlation coefficient of 500 hPa temperature in the North Hemisphere during 1-15 May 2013 (a) anomaly correlation coefficient, (b) significance level

以上的温湿及其倾向廓线的分析结果均显示:C-P边界层方案避免了物理过程调用前后的插值，减小了插值误差和插值带来的振荡，较原来的Lorenz跳点调用更为合理。需要说明的是, 由于水平风在半层，无需插值，因此, 控制试验中水平风及其倾向的廓线分布未发现异常。这也说明边界层过程调用前后的插值极有可能是造成温湿廓线锯齿状振荡的重要原因，而采用C-P边界层方案避免耦合时的插值是改进GRAPES全球模式边界层过程的有效手段。连续试验结果也表明, C-P边界层方案应用在一定程度上提高了模式预报的综合性能。

将Lorenz跳点计算的物理过程温湿变量通过插值与C-P跳点的动力框架耦合, 会造成边界层物理过程温湿廓线的异常, 下面分析其原因。式(1)中忽略反梯度输送γ，并假设交换系数K为常数，问题简化为一维热传导方程的数值求解，方程右端为ψ的二阶导数。给出定解问题：

(7)

式(7)中，T表示温度，边条件为T(z=0)=T(z=H)=0，初条件设为T(z, 0)=bz(H-z)/H²。这是一个两端温度为0，无热源的热传导问题。取H=1000，a=1，b=1。

利用高时空分辨的时间向前显式积分数值结果作为可参考的真值解。设计以下几组数值求解方案进行对比：

①时间向前显示积分，整层T用线性插值到半层，求出半层倾向后，用线性插值将半层的倾向插回到整层，更新T求解。

②隐式积分，整层直接求解(该方案与C-P边界层方案与框架耦合流程一致)。

③利用线性插值后，在半层隐式积分，求出半层倾向后，用线性插值再插回到整层，更新T求解(这种方案和控制试验中GRAPES全球模式中物理过程和动力框架耦合流程一致)。

在①和③方案中，数值求解的第2类边条件设为。各种方案在2 h的倾向值和积分数值解见图 7。图 7中黑实线是高分辨积分结果，可看作真解。由图 7可以看出，用隐式积分方案直接在整层上计算其结果和显示方案在整层计算的结果非常接近。在半层上用线性插值的显示积分结果稍微有一些误差，但其趋势与整层计算结果相近。而采用隐式积分，同时利用线性插值在半层上的计算结果出现了大幅度的振荡。这种振荡与GRAPES_GFS边界层中出现的振荡非常相似。由此推测通过线性插值进行边界层物理过程和动力框架之间的耦合极有可能是造成振荡的原因。

图7

图7 不同方案积分2 h后的结果(a)整层上的瞬时温度倾向，(b)温度 Fig.7 Forecast results at 2 h integration (a) temperature tendeny at the full level, (b) temperature

通过对一维热传导方程的数值分析发现，两次插值过程在隐式求解方程时导致了温度的锯齿状振荡，而直接在整层求解避免了插值误差，与高分辨积分结果几乎一致，与插值的耦合方案比较，其精度最高。可见在物理过程和动力框架耦合过程中，避免插值对提高模式的精度相当重要。

通过对边界层热量交换系数的计算和下边界条件处热通量常值通量层的考虑，进行了C-P边界层方案的研发，并将其应用于GRAPES全球模式。从个例试验和连续试验的分析可知：

1)由于GRAPES全球模式中采用温度和湿度的插值将边界层过程和动力框架耦合，而插值方案的二阶不光滑极有可能导致边界层方程中含有的二阶导数的隐式积分不协调，从而造成边界层温湿倾向和边界层中温湿廓线锯齿状振荡。

2) GRAPES全球模式中设计了C-P边界层方案，避免了插值对边界层方程中二次导数求解的影响，有效消除了模式预报的边界层中温度和湿度不合理的锯齿状振荡，使模式中边界层过程的描述趋于合理。

3)从结果分析看，C-P边界层方案减小了模式的预报误差，提高了模式的整体效果。这主要是由于C-P边界层方案避免了插值误差，使耦合精度提高。另一方面，由于模式的垂直方向是不等距分层，因此, 在C-P跳点上求解温湿的一阶导数相对于Lorenz跳点精度更高。尽管二阶导数产生了距离误差，但是相比较而言，一阶导数是大值，在差分计算中更重要，这也是在C-P跳点上计算边界层温湿倾向的优势。

C-P边界层方案的设计原理中提到，在稳定边界层中的局地方案中，仍采用插值计算局地R_i数，因此, 插值不可能完全避免。但边界层过程中对于大尺度反馈起决定性作用的主要部分是不稳定边界层中的湍流交换，而C-P边界层方案中不稳定边界层的交换系数直接求取，无需插值；且C-P边界层方案避免了边界层过程和动力框架耦合时整个垂直层的两次插值。因此, 从理论分析，GRAPES全球模式中使用C-P边界层方案比原来控制试验中所用方案优越。但采用C-P跳点由于水平风和温度不在同一层，在计算边界层高度和局地R_i数等参数时较为复杂，需要考虑相关变量的插值。C-P边界层方案是针对我国GRAPES模式的框架特点开发，因此, 对于国际通用的研究模式如WRF其应用价值有局限性。从综合效果看，C-P边界层方案对GRAPES模式只是稍有改进，毕竟边界层过程只是模式中的一个部分，也期待其他物理过程和动力框架的耦合方案得到改善，以进一步提高模式预报的精度和综合性能。