Possible units error in getVerticalMassFlux

[sdeastham]

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!

[mathomp4]

I'm going to assign @wmputman and hope he sees this. It's beyond me!

[tclune]

Yes - ultimately this is an issue that only @wmputman can answer. If we don't hear back in a timely manner, we can at least have someone repeat this analysis as a check.

[sdeastham]

@tclune - think that @wmputman is not able to look at this - any thoughts on next steps?

[wmputman]

I have to look at it, sorry for the delay

[wmputman]

This was taken long ago from the old LatLon core. This code is indeed wrong. You just need to remove the fac in the calculation of conv:

MFX = Pa m2 / s MFY = Pa m2 / s

conv & pit will maintain units of Pa * m2 / s

(conv - bkpit) / (grav area) = ( Pa m+2 s-1 ) / ( m s-2 m+2 ) = ( N s-1 ) / ( m+3 s-2 ) = N s+1 m-3 = kg m-2 s-1

If you have a bad MFZ, you should be able to correct it by multiplying by gravity*dt (dt=450s)

[mathomp4]

@wmputman Queries:

1. Does this affect GEOS or is this in a code path we don't use?
2. Do you have a possible fix for the code? Or I guess @sdeastham could propose one.

[wmputman]

It only changes the exports for GEOS. I have committed a code fix:

https://github.com/GEOS-ESM/FVdycoreCubed_GridComp/commit/99b2e9a123525bd18fa366d778f852530101bb64#diff-3dd8e979b916d8034f7682465fb7f31ce0021dda0299b0314edaa496387a21df

[sdeastham]

Thanks @wmputman for the review and fix! @lizziel - suggest we cherry pick at least this update into the GCHP fork, since it is trivial but will make the diagnostic much more useful?