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:

FVdycoreCubed_GridComp/DynCore_GridCompMod.F90

Line 4347 in b331e2a

call getVerticalMassFlux(mfxxyz, mfyxyz, mfzxyz, dt)

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:

FVdycoreCubed_GridComp/DynCore_GridCompMod.F90

Line 998 in b331e2a

UNITS = 'Pa m+2 s-1', &

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:

FVdycoreCubed_GridComp/DynCore_GridCompMod.F90

Line 1032 in b331e2a

UNITS = 'kg m-2 s-1', &

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

FVdycoreCubed_GridComp/FV_StateMod.F90

Lines 3178 to 3179 in b331e2a

	conv(i,j,k) = ( xfx(i,j,k) - xfx(i+1,j,k) + &
	yfx(i,j,k) - yfx(i,j+1,k) ) * fac

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:

FVdycoreCubed_GridComp/FV_StateMod.F90

Line 3170 in b331e2a

fac = 1.0/(dt*MAPL_GRAV)

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:

FVdycoreCubed_GridComp/FV_StateMod.F90

Line 3204 in b331e2a

mfz(i,j,k) = ( conv(i,j,k-1) - FV_Atm(1)%bk(k)*pit(i,j) )/(MAPL_GRAV*fv_atm(1)%gridstruct%area(i,j)) ! Kg/m^2/s

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 - bk*pit) / (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:

99b2e9a#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?