Estimate diffuse and direct component from global irradiance

I am looking to separate the diffuse and direct component of global irradiance and found the Erbs model to do this in pvlib (see pvlib.irradiance.erbs) however, I am getting very strange results. I would expect the Direct Normal Irradiance (DNI) to be lower than the Global Horizontal Irradiance (GHI); or am I missing something? Values of GHI are not above 800 W m^2 for these days.

EDIT: As per Cliff H advice, I have limited the solar zenith to less than 85 arc degrees; the results have improved however, there are large spikes in DNI values that do not appear very reasonable, e.g. start of 07-16.

DNI > GHI is common at low solar elevation. GHI decreases much faster than DNI as solar elevation decreases. For example, think of a clear day with the sun right near the horizon. DNI will large because it's measured on a plane normal to the sun vector, but GHI will be near zero.

The values of DNI that are much greater than 1000 W/m2 are likely at very high zenith, since the Erbs model basically divides by cos(zenith). In practice, I limit using decomposition models like Erbs to zenith<85 degrees, to avoid the non-physical results.