Yes, indeed. Great catch! The issue here is that the time-averaged output NetCDF files (when AVERAGES or DIAGNOSTICS_* are activated) are delayed in its creation until the fields are accumulated and ready to write out if the parameters NDEFAVG or NDEFDIA are greater than zero (multi-file option).
It is tricky, but the solution is very simple. The AVG(ng)%load or DIA(ng)%load are reset to zero for such cases in read_phypar.F. By default, their values are initialized to one during initialization in load_s1d. It is weird that the compiler didn't pick-up the out-of-bounds errors during my original testing.
