WRF vertical integrated moisture flux using NCL
I am going to explain my answer based on how the authors from this paper calculate this Vertically integrated moisture flux convergence as a predictor of thunderstorms. This is the equation that is defined in their paper

$$ VIMFC = -\frac{1}{g} \int_{700\,hPa}^{1000\,hPa}\ (\frac{\partial u q}{\partial x}+\frac{\partial v q}{\partial y}) * dp $$

For a reference for the vertically integrated moisture flux please look at the equation in the reference A Hydrological Definition of Indian Monsoon Onset and Withdrawal

$$ VIMFT = -\frac{1}{g} \int_{700\,hPa}^{1000\,hPa} (q * u + q * v) * dp$$

So I will provide an algorithm for a numerical implementation of this equation and opinions can vary what exactly constitutes a "atmospheric layer" in NWP. 

In this equation q is the specific humidity, u and v are the zonal and meridional velocities, p is the pressure,  g is the acceleration due to gravity. 

Moisture flux convergence or divergence is 

$$ (\frac{\partial u q}{\partial x}+\frac{\partial v q}{\partial y})$$

Moisture flux is 

$$ q * u + q * v  $$

You then need to sum the horizontal moisture flux convergence(or in your case moisture flux) over the following layers
1) 1000 - 925 hPa
2) 925-850 hPa
3) 850-700 hPa
So the definition of what a "layer" is subjective and you can include three pressure surfaces in it if you think a deep layer is required as opposed to a shallow layer. Numerically a layer can be thought as the average values of two pressure surfaces. As an example in your case you take the average of the velocities at the 1000 hPa and 925 hPa surfaces and the average of the specific humidity at the 1000 hPa and 925 hPa.

Once you calculate the horizontal moisture flux(or if you are calculating MFC then you need to use  finite centered differences and backward and forward differences at the grid boundaries and that procedure is explained in this ESSE answer - Moisture flux convergence numerical procedure and here Finite differences for calculating MFC) over each layer you multiply that by dp and then sum up the contributions. 

UPDATE

In response to the questions in the comments 

$"1000 - 925"$

is not to be interpreted as a subtraction but as an average.

What needs to be done is to take an average of two pressure surfaces in the following manner

$ q_{layer} = (q_{1000} + q_{925})/2$

$ u_{layer} = (u_{1000} + u_{925})/2$

$ v_{layer} = (v_{1000} + v_{925})/2$

Once you do this for every layer multiply by the corresponding value of dp and in this case 

$ dp = 1000 -925 $

This is a subtraction. Then you plug the individual values into the sum of individual contributions for each layer and then sum the layers

$ (q_{layer}* u_{layer} + q_{layer}* v_{layer}) * dp$

The units of VIMFC is typically 

$$ 10^− 5 \, kg \, m^− 2 \,s^− 1 $$