Potential Cancellation Issue in deln_flux

[oelbert]

The d2 calculation near the end of deln_flux has discrepant results between 64- and 32-bit caulculations:

GFDL_atmos_cubed_sphere/model/tp_core.F90

Line 1316 in 9ca1468

d2(i,j) = (fx2(i,j)-fx2(i+1,j)+fy2(i,j)-fy2(i,j+1))*gridstruct%rarea(i,j)

The GEOS team has resolved this by enforcing the order-of-operations in the flux addition/subtractions:
d2(i,j) = ((fx2(i,j)-fx2(i+1,j))+(fy2(i,j)-fy2(i,j+1)))*gridstruct%rarea(i,j)

This could be a good improvement to 32-bit runs

[lharris4]

Hi, @oelbert. Thank you for identifying these changes from GMAO. Indeed this better representation of order-of-operations when computing flux divergences, which are done at a number of places in the solver, could go a long way to avoiding some of the cancellations that can occur.

They also define the Flux Capacitors mfx,mfy and accumulated Courant numbers cx and cy as double-precision in all runs, and a double-precision version of delp. I definitely think these are worth trying out.