Possible units error in getVerticalMassFlux

I have been comparing the estimated vertical mass fluxes generated by FV3 with those coming from the older FV code embedded in GEOS-Chem, and have found significant disagreement. However, I am wondering if this is because of a units issue and would appreciate any insights that the developers might have.

Specifically, in FV_StateMod, the vertical mass flux mfzxyz is calculated by calling fv_getVerticalMassFlux with mfxxyz and mfyxyz as inputs:
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/DynCore_GridCompMod.F90#L4347

Immediately beforehand, mfxxyz and mfyxyz are used to fill the MX and MY exports, which are stated to be in Pa m+2 s-1:
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/DynCore_GridCompMod.F90#L998

Similarly, mfxxyz is used immediately after the call to fv_getVerticalMassFlux to fill the export MFZ, which has declared units of kg m-2 s-1:
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/DynCore_GridCompMod.F90#L1032

However, as far as I can tell, the operations in fv_getVerticalMassFlux will not convert a quantity with units Pa m+2 s-1 to a quantity with units kg m-2 s-1. Briefly, it appears that the routine in question first calculates conv, which must have units of Pa m+2 s-1 multiplied by the units of fac (noting that xfx = mfx):
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/FV_StateMod.F90#L3178-L3179

fac is equal to 1.0/(dt*MAPL_GRAV), and must therefore have units of s-1 * s+2 * m-1 = s+1 m-1:
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/FV_StateMod.F90#L3170

That implies that conv has units of Pa m+2 s-1 * s+1 m-1 = Pa m. The only remaining operation in fv_getVerticalMassFlux which should affect the units of the answer are on line 3204, when mfz is calculated. Here, b(k) * pit is subtract from conv, and then the result divided by MAPL_GRAV*area:
https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/blob/b331e2ae0339bde84d63eacd911c36befe74d735/FV_StateMod.F90#L3204

pit is just an accumulation of conv so should have the same units, which I believe are still Pa m. Dividing by MAPL_GRAV*area should then have the effect of changing the units to Pa m * s+2 m-1 * m-2 = Pa s+2 m-2. Converting from pressure to mass still gives Pa = N m-2 = kg m s-2 m-2 = kg s-2 m-1, which means the eventual units end up being kg m-3. That would imply that a factor with units equal to velocity is missing from the calculation.

I haven't been able to find a fundamental error in my calculation, but may well be missing something obvious. Nonetheless, I would appreciate any insight that anyone can provide. My suspicion (but it is that at most) is that the application of fac is incorrect - removing that factor at least causes the units to work out correctly. It also makes logical sense, since inclusion of fac results in a duplicate division by MAPL_GRAV, and incurs a division by dt when the units of mfx are already meant to be a tendency. Alternatively, it may be that the units of mfx and mfy are incorrectly listed; but again, any information anyone can provide would be very welcome!