【数值模型后处理系列】通风系数计算及垂直层插值

一、通风系数

1.1 通风系数简介

通风系数（Ventilation Coefficient，VC）可以用来表征扩散条件，其计算公式如下（参考U S Iyer and P Ernest Raj的文章）：
在这里插入图片描述
其中mixing depth选用WRF输出的边界层高度（PBLH），mean wind speed近似用边界层顶的风速与地面风速做平均（当然也可多选几层）。

1.2 Python代码实现VC的计算

计算VC的示例代码：

from netCDF4._netCDF4 import Dataset
from wrf import getvar, ALL_TIMES, interplevelfiles = 'wrfout文件列表'
wrfout = [Dataset(i) for i in files]pblh = getvar(wrfout, 'PBLH', timeidx=ALL_TIMES, method='cat')
height = getvar(wrfout, "height_agl", timeidx=ALL_TIMES, method='cat')
wind = get_uvmet_wspd_wdir(wrfout, timeidx=ALL_TIMES, method='cat')
wind_pblh = interplevel(wind[0], height, pblh)
vc = (wind_pblh + wind[0, :, 0, :, :]) / 2 * pblh  # ventilation coefficient(m2/s)

这样计算出来的VC会有许多网格点是空值，可以在interplevel中指定missing参数，不过只能指定一个数值，由于各个网格的VC数值并不知道也不相等，因此，可以插值后再进行缺失值补空，这部分有许多方法可以实现，读者可自行研究。

1.3 结果解读

计算出来的VC空间分布图如下（川渝地区）。由图可知，该时刻，川渝地区通风系数较小，扩散条件相对较差。

在这里插入图片描述

二、垂直层差值

2.1 插值函数简介

上面计算VC需将风速插值到PBL高度层，wrf-python中提供了interplevel函数可满足这一插值需求，

wrf.interplevel(field3d, vert, desiredlev, missing=<MagicMock name='mock().item()' id='140675643013392'>, squeeze=True, meta=True)

参数解释：

field3d (xarray.DataArray或numpy.ndarray) – 要插值的三维字段，最右边的维度为 nz × ny × nx；
vert (xarray.DataArray或numpy.ndarray) – 垂直坐标的三维数组，通常是压力或高度。该数组必须具有与field3d相同的维度；
desiredlev ( float、一维序列或numpy.ndarray) – 所需的垂直水平。这可以是单个值（例如 500）、值序列（例如 [1000, 850, 700, 500, 250]）或多维数组，其中右侧二维 (ny x nx) 必须匹配field3d，并且任何最左边的尺寸与 field3d.shape[:-3] 匹配（例如行星边界层）。必须与vert参数采用相同的单位；
missing ( float) – 用于输出的填充值。默认为wrf.default_fill(numpy.float64);
squeeze（bool，可选） – 设置为 False 以防止大小为 1 的维度自动从输出形状中删除。默认为 True；
meta ( bool) – 设置为 False 以禁用元数据并返回 numpy.ndarray而不是 xarray.DataArray.默认为 True。

函数返回：插值变量。如果启用 xarray 并且meta参数为 True，则结果将是一个 xarray.DataArray对象。否则，结果将是一个numpy.ndarray没有元数据的对象。

返回类型：xarray.DataArray or numpy.ndarray

2.2 使用示例

示例1：将位势高度插值至 500 hPa

from netCDF4 import Dataset
from wrf import getvar, interplevelwrfin = Dataset("wrfout_d02_2010-06-13_21:00:00")p = getvar(wrfin, "pressure")
ht = getvar(wrfin, "z", units="dm")ht_500 = interplevel(ht, p, 500.0)

示例2：将相对湿度插值到边界层高度

from netCDF4 import Dataset
from wrf import getvar, interplevelwrfin = Dataset("wrfout_d02_2010-06-13_21:00:00")rh = getvar(wrfin, "rh")
height = getvar(wrfin, "height_agl")
pblh = getvar(wrfin, "PBLH")rh_pblh = interplevel(rh, height, pblh)