问如何计算Python中水汽通量辐合的垂直积分？

EN

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
import pandas as pd
import cartopy.feature as cf
from netCDF4 import Dataset 
import metpy.calc as mpcalc

ds1 = xr.open_dataset('Downloads/uwnd.mon.mean.nc')
ds2 = xr.open_dataset('Downloads/vwnd.mon.mean.nc')
ds3 = xr.open_dataset('Downloads/shum.mon.mean.nc')
lons = ds1['lon'][:]  
lats = ds1['lat'][:] 
levs = ds1['level'][2]

# --- select all JJAS months
uJJAS_ = ds1['uwnd'].sel(time=np.in1d(ds1['time.month'], [6, 7, 8,9]))
uJJAS_
uJJAS=uJJAS_[:,2,:,:]
uJJAS
# Select JJAS 850 hPa from 1981-2010    
u_1=uJJAS[132:252,:,:]
u_1
# Prepare the climatology
u_1_mean= np.mean(u_1,axis=0)
u_1_mean

# --- For v
vJJAS_ = ds2['vwnd'].sel(time=np.in1d(ds2['time.month'], [6, 7, 8,9]))
vJJAS_
vJJAS=vJJAS_[:,2,:,:]
vJJAS    
v_1=vJJAS[132:252,:,:]
v_1
v_1_mean= np.mean(v_1,axis=0)
v_1_mean

# --- For specific humidity
shJJAS_ = ds3['shum'].sel(time=np.in1d(ds3['time.month'], [6, 7, 8,9]))
shJJAS_
shJJAS=shJJAS_[:,2,:,:]
shJJAS    
sh_1=shJJAS[132:252,:,:]
sh_1
sh_1_mean= np.mean(sh_1,axis=0)
sh_1_mean

uq=(u_1_mean)*(sh_1_mean)
vq=(v_1_mean)*(sh_1_mean)

# Compute dx and dy spacing for use in divergence calculation
dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)

# Calculate moisture flux divergence
HMC_LE = (np.array(mpcalc.divergence(uq, vq, dx=dx, dy=dy)))

lat=lats.to_numpy()
lon=lons.to_numpy()
u_=uq.to_numpy()
v_=vq.to_numpy()
ax2 = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
clevs = [-0.06,-0.05,-0.04,-0.03,-0.02,-0.01, 0,0.01,0.02,0.03,0.04,0.05,0.06]
vp_fill = ax2.contourf(lons,lats,HMC_LE*1000,clevs,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r,
extend='both')
plt.colorbar(vp_fill, orientation='vertical')
q2=ax2.quiver(lon, lat, u_, v_,width=0.003,headwidth=5, scale_units='xy',scale=25,transform=ccrs.PlateCarree())
ax2.coastlines(alpha=0.8)
ax2.set_xticks([60,70,80,90,100], crs=ccrs.PlateCarree())
ax2.set_yticks([-10,0,10,20,30,40], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,number_format='.0f')
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)
ax2.set_extent([57,101,-10,45])
ax2.add_feature(cf.BORDERS,edgecolor='black',alpha=0.8)
plt.title('Moisture flux divergence at 850 hPa ', fontsize=16)
plt.show()