用Python计算湿位涡错误

逆光--

通过公众号（好奇心Log）的推文改写的计算湿位涡剖面的程序，与文献使用同样的资料求的是2016年8月21日20时（b）沿106°E的MPV1 经向剖面图,差异也太大了，是程序计算哪里有问题吗？求大神指教！

import cfgrib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import metpy.calc as mpcalc
from metpy.interpolate import cross_section
from metpy.units import units
from metpy.constants import g
import scipy.ndimage as ndimage

data = cfgrib.open_dataset('......\gdas1.fnl0p25.2016082112.f00.grib2',filter_by_keys={'typeOfLevel': 'isobaricInhPa'}).metpy.parse_cf()
data

heights = data['gh']
ug, vg = mpcalc.geostrophic_wind(data['gh'])

vert_abs_vort = mpcalc.absolute_vorticity(ug, vg)
vert_abs_vort

temperature, pressure, relative_humidity = xr.broadcast(data['t'],data['isobaricInhPa'],
                                                        data['r'])
#露点温度
dewpoint = mpcalc.dewpoint_from_relative_humidity(temperature,relative_humidity)
dewpoint

# 相当位温
theta_e = mpcalc.equivalent_potential_temperature(pressure, temperature, dewpoint)
theta_e

data['theta_e'] = xr.DataArray(theta_e.metpy.magnitude,
                               coords=heights.coords,
                               dims=heights.dims,
                               attrs={'units': 'K'})
data['u_g'] = xr.DataArray(ug.metpy.magnitude,
                           coords=heights.coords,
                           dims=heights.dims,
                           attrs={'units': 'm/s'})
data['v_g'] = xr.DataArray(vg.metpy.magnitude,
                           coords=heights.coords,
                           dims=heights.dims,
                           attrs={'units': 'm/s'})
data['relative_humidity'] = xr.DataArray(relative_humidity.metpy.magnitude,
                           coords=heights.coords,
                           dims=heights.dims,
                           attrs={'units': 'percent'})
data

dtheta_e_dp, dtheta_e_dy, dtheta_e_dx = (var.metpy.unit_array for var in mpcalc.gradient(data['theta_e'], axes=('vertical', 'y', 'x')))
dug_dp = mpcalc.first_derivative(data['u_g'], axis='vertical')
dvg_dp = mpcalc.first_derivative(data['v_g'], axis='vertical')
dz_dp = mpcalc.first_derivative(data['gh'], axis='vertical')
#mpv = g * (1 / dz_dp) * (-dvg_dp * dtheta_e_dx + dug_dp * dtheta_e_dy + vert_abs_vort * dtheta_e_dp)
#mpv
mpv1 = g * (1 / dz_dp) * (vert_abs_vort * dtheta_e_dp)
mpv2 = g * (1 / dz_dp) * (-dvg_dp * dtheta_e_dx + dug_dp * dtheta_e_dy)
mpv1

data['mpv1'] = xr.DataArray(mpv1*10**4,
                           coords=heights.coords,
                           dims=heights.dims,
                           attrs={'units': 'microkelvin/s^3'})
data['mpv2'] = xr.DataArray(mpv2*10**4,
                           coords=heights.coords,
                           dims=heights.dims,
                           attrs={'units': 'microkelvin/s^3'})
data

# Cross section parameters
start = (35.00,106.00)
end = (45.00,106.00)

cross = cross_section(data,start,end)
cross

fig = plt.figure(1, figsize=(14, 6))
ax = plt.axes()

#mpv1_contour = ax.contourf(cross['latitude'], cross['level'], cross['mpv1'], extend='both',
#                        levels=np.arange(-1, 0.2, 0.2), cmap='bwr')

#mpv1_colorbar = fig.colorbar(mpv1_contour)

mpv1_contour = ax.contour(cross['latitude'], cross['isobaricInhPa'], cross['mpv1'],
                           levels=np.arange(-3.0, 2.0, 0.5), colors='k', linewidths=1,
                           linestyles='-', alpha=1, zorder=2)
mpv1_labels = ax.clabel(mpv1_contour,inline=True,fmt='%.0f',fontsize=10,colors='tab:blue')

##mpv2_contour = ax.contour(cross['latitude'], cross['level'], cross['mpv2'],
#                           levels=np.arange(0, 1, 0.2), colors='k', linewidths=1,
#                           linestyles='-', alpha=0.5, zorder=2)
##mpv2_labels = ax.clabel(mpv2_contour,inline=True,fmt='%.0f',fontsize=10,colors='k')

ax.set_yscale('symlog')
ax.set_yticklabels(np.arange(1000, 100, -100))
ax.set_ylim(1000,200)
ax.set_yticks(np.arange(1000, 100, -100))

plt.show()

文献图b

python画的结果

雪域小丹

是不是积分问题？后来解决成功了吗？求共享！

逆光--

雪域小丹 发表于 2023-2-9 12:03
是不是积分问题？后来解决成功了吗？求共享！

还没有解决

ququ

楼主解决了吗？求共享

tulalang

楼主你程序里只绘制出了mpv1呢，会不会是应该绘制mpv1+mpv2，没数据也很难帮你验证

左.

为什么这里面的湿位涡计算的符号和公式里的相反啊，求教

逆光--

tulalang 发表于 2023-4-7 08:37
楼主你程序里只绘制出了mpv1呢，会不会是应该绘制mpv1+mpv2，没数据也很难帮你验证

用的是NCEP再分析资料(水平分辨率0.25°×0.25°)2016年8月21日20时的数据。只绘制mpv1是因为跟着文献绘制的。（破涕为笑.jpg）

逆光--

左. 发表于 2023-4-13 12:52
为什么这里面的湿位涡计算的符号和公式里的相反啊，求教

dz_dp是负值