Reg, I think I now have a more accurate working program.
I am computing the E and H fields at only 20m, which
produces dramatically different results. Using
the Poynting Vector, S = 1/2*Re(E X H*), where
"X" is the vector cross product (H* is the complex conjugate
of H). Integrating |S| over a hemispherical region, of radius
20 m, shows a radiation efficiency of 80.3%, and a radiation
resistance of 27.7 ohms.

I have moved the test frequency to 8.1 MHz where the input
impedance of the antenna (99 x 10 m radials) is:
34.523 +j0.18. This implies a radial input impedance of
6.823 + j 0.18 ohms. These results appear to be much closer
to your program. Of course, now I am starting to wonder
if I should not redo the computation at 10 m to see how
significant the ground wave losses are. Closer than 10 m
is probably not practical since I do not think NEC can
compute the near fields in cylindrical coordinates.

The results also verify your comments concerning the
contribution of the ground wave to the total radiated
power. The "Sky wave" radiated power represents
only about 35% of the total input power.

My new Excel spread sheet contains over 4,000 active
cells, but is much easier to use than the previous method.
Using rotational symmetry, and other methods, the
NEC run time has been dramatically reduced.
I have also used almost lossless wire of conductivity
1E12 S/m. Perfect wire crashes NEC when using
the "Numerical Green's Function" -- which helps speed
up calculations.

To verify program accuracy I have computed the radiation
efficiency of an ideal 9m monopole over a perfectly conducting
ground. The results are accurate to within less than 0.25%.

Obviously my previous results are no longer valid, So will
have to re-calculate all the test antennas.

Frank