在 Python 中从 WRF 输出计算 3D Frontogenesis

Question

我目前正在尝试运行 python 脚本来计算我的 WRF 输出数据上的 3 维锋生，以执行横截面分析。我已经有一个可用的 Petterson 公式的二维版本来进行比较，但是我的三维版本只捕捉了二维版本的非常小的方面。

我用来计算 3D 锋生的公式是：

其中梯度项为：

这是我的二维码 (925 hPa) 生成的示例图像：

这是一个三维代码，然后向下插值到相同的二维 (925 hPa) 表面：

一些读数至少出现在图像上的相同地理区域这一事实向我表明，我至少在正确的轨道上，并且我可能在某个地方犯了一个小错误。从我所见，我的 python 代码看起来是正确的，至少根据我对 np.gradient 函数如何工作的理解，这里是计算前沿生成的代码：

# Fetch the fields we need
p = getvar(ncFile, "pressure") * 100 # Convert to Pa
z = getvar(ncFile, "z")
ua = getvar(ncFile, "ua")
va = getvar(ncFile, "va")           
theta = getvar(ncFile, "theta")
omega = getvar(ncFile, "omega")

dz = np.gradient(z, axis=0)
dp = np.gradient(p, axis=0)

theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (np.gradient(theta, axis=0) / dz)**2)

zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1)))
meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1)))
vertical_gradient = (-1 * (np.gradient(theta, axis=0) / dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1)))

F3D = 1.08e9 * (1 / theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient)
return F3D

作为参考，dx 和 dy 项是直接从 NetCDF 文件本身获取的（它有属性 DX 和 DY，都定义为 4000m）

我正在使用 wrf-python 库通过 getvar 导入我的数据，它正在导入 netCDF 文件。作为参考，netCDF 文件使用类似于标准 numpy 数组的数组结构： [自下而上、南北、东西]

所以轴参数的顺序应该是正确的（z = 0，y = 1，x = 2）。

我的一位指导老师认为梯度的边缘可能会导致计算内部出现一些问题，但我的理解是每个点都是独立于边缘计算的，因此应该没有区别，但这是我不是 100% 确定的。

有谁知道为什么计算会产生如上图所示的错误结果？

Answer 1

在多打了几下之后，我能够解决这个问题。去展示以确保考虑您的完整单元分析！

位温梯度中的项是考虑纵轴为z，而垂直梯度中关于欧米茄的项是考虑p。由于欧米茄是压力的垂直单位，因此我对潜在温度梯度的术语选择不正确。将导数从 z 换成 p 解决了问题的前半部分。

其次，在计算垂直方向的导数时，您需要考虑到 numpy.gradient 假设您已经经过两个数据点，因此它错误地将结果除以二，所以我创建了一个函数包装器来处理分析的内点和外点：

def calc_center_difference(A, ax):
    gradient = np.gradient(A, axis=ax) 
    gradient *= 2.0
    if ax == 0:
        gradient[0,:,:] /= 2.0
        gradient[-1,:,:] /= 2.0   
    elif ax == 1:
        gradient[:,0,:] /= 2.0
        gradient[:,-1,:] /= 2.0
    elif ax == 2:
        gradient[:,:,0] /= 2.0
        gradient[:,:,-1] /= 2.0
    else:
        return ValueError("Invalid axis passed to calc_center_difference")
    return gradient

之后，我将考虑 dz 或 dp 的衍生术语的实例与包装函数实例交换，以确保这些术语被正确计算，瞧！一切正常！

这是我之前的函数的更正形式，可以在 python 中计算完整的 3D 锋生：

# Input netcdf
# - [bottom_top, north_south, west_east]
def three_d_fronto(ncFile):
    # Fetch the fields we need
    p = to_np(getvar(ncFile, "pressure") * 100)
    z = to_np(getvar(ncFile, "z"))
    ua = to_np(getvar(ncFile, "ua"))
    va = to_np(getvar(ncFile, "va"))    
    theta = to_np(getvar(ncFile, "theta"))
    omega = to_np(getvar(ncFile, "omega"))

    dz = calc_center_difference(z, 0)
    dp = calc_center_difference(p, 0)

    theta_gradient = np.sqrt((np.gradient(theta, dx, axis=2))**2 + (np.gradient(theta, dy, axis=1))**2 + (calc_center_difference(theta, 0) / dp)**2)
    zonal_gradient = (-1 * np.gradient(theta, dx, axis=2)) * ((np.gradient(ua, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dx, axis=2) * np.gradient(theta, dy, axis=1)))
    meridional_gradient = (-1 * np.gradient(theta, dy, axis=1)) * ((np.gradient(ua, dy, axis=1) * np.gradient(theta, dx, axis=2)) + (np.gradient(va, dy, axis=1) * np.gradient(theta, dy, axis=1)))
    vertical_gradient = (-1 * (calc_center_difference(theta, 0) / dp)) * ((np.gradient(omega, dx, axis=2) * np.gradient(theta, dx, axis=2)) + (np.gradient(omega, dy, axis=1) * np.gradient(theta, dy, axis=1)))

    F3D = 1.08e9 * (1 / theta_gradient) * (zonal_gradient + meridional_gradient + vertical_gradient)
    return F3D

提醒一下，您需要安装 wrf-python 库才能使用此功能。