How metpy module calculate Moisture flux divergence in python

I am calculating moisture flux divergence using metpy module in python. I calculated it like `

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import pandas as pd
import metpy.calc as mpcalc
import matplotlib.patches as mpatches
file1=xr.open_mfdataset('/2018,19_1000to500/*.nc')
file=file1.sel(longitude=slice(75.75, 77.25),latitude=slice(10.5,9.5))
q1=(file.q)*1000
u1=file.u
v1=file.v
lat=file.latitude
lon=file.longitude
for k in range (1):  
    level=['1000'] #,,'925','850','700','600','500'
    
    q1d=q1.resample(time='1D').mean().sel(level=int(level[k]))
    u1d=u1.resample(time='1D').mean().sel(level=int(level[k]))
    v1d=v1.resample(time='1D').mean().sel(level=int(level[k]))
    sdate=['2018-08-08','2018-10-06','2018-11-10','2019-06-06','2019-08-06','2019-10-24'] #,
    edate=['2018-08-17','2018-10-15','2018-11-19','2019-06-13','2019-08-14','2019-10-27']#,,,
    for i in range(6):
        q1dm=q1d.sel(time=slice(sdate[i],edate[i])).mean('time')
        u1dm=u1d.sel(time=slice(sdate[i],edate[i])).mean('time')
        v1dm=v1d.sel(time=slice(sdate[i],edate[i])).mean('time')
        qu=q1dm*u1dm
        qv=q1dm*v1dm
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
 HMD = (np.array(mpcalc.divergence(qu, qv, dx=dx, dy=dy))) 
 print (HMD)
        

` but I knew that metpy calculate divergence as the sum of du/dx and dv/dx. So my question is if I put qu,qv instead of u and v in the code how metpy calculate the moisture divergence. Will it calculate it as the sum of dqu/dx and dqv/dy.

I manually calculated the sum of dqu/dx and dqv/dy and compared it with the values provided by metpy and there is very very small difference in it(the difference starts in the 3rd or 4th decimal place). I am asking this because the documentation of metpy moisture divergence is not available anywhere , So I need to know how metpy compute it what is the exact equation it uses for the computation.

MetPy is calculating divergence as du/dx + dv/dy. So if you pass in u=qu and v=qv, then the return from divergence will be dqu/dx + dqv/dy.

The small differences you note may be due to how MetPy is calculating du/dx and dv/dy, which accounts for issues around calculating derivatives on the spherical surface of the Earth/projected grid and handles non-uniform grid spacing in its finite difference formulation in approximating the derivatives.