How to calculate Vertical Integration of moisture flux convergence in Python?

Question

I'm new to python, I calculated the Moisture Flux Divergence [MFD] (i.e. {d(qu)/dx}+{d(qv)/dy} ,q is specific humidity and u,v are zonal and meridional wind components) over India at pressure level 850 hPa using metpy.calc function HMC_LE = (np.array(mpcalc.divergence(uq, vq, dx=dx, dy=dy))). I have followed the steps described here: How to calculate moisture flux divergence in python. My code -

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()

The resulting plot is shown here - I have values for 1000 hPa, 925 hPa, 850 hPa, 700 hPa, 600 hPa, 500 hPa, 400 hPa, 300 hPa. Now I want to integrate MFD values vertically [from 1000 hPa to 300 hPa pressure level]. Please tell how to do this integration? (INTEGRAL from 300 hPa to 1000 hPa of ({d(qu)/dx}+{d(qv)/dy}) dP.)

Thank you!

Answer 1

If you assemble all the levels into a single array, you should be able to integrate using numpy.trapz or some of the more complicated integration methods in scipy.integrate.