I implemented the standard numerical algorithm for calculating the Berry curvature in MATLAB. For a given 2D system, I can visualize the Berry curvature over parameter space. If I sum the Berry curvature at all points in a region of parameter space, I am effectively taking the Berry curvature area integral inside this region. If this region is the entire Brillouin zone, I get the first Chern number.
Due to Stokes' theorem, this area integral is nothing but the gauge-invariant Berry phase due to adiabatic traversal along a parameter space loop bordering this region. In other words, in the absence of gauge singularities, the Berry curvature $F$ area integral is equal to the Berry connection $A$ line integral: $$ \oint A = \iint F. $$
For my question, consider a graphene-like system where the energy dispersion is degenerate at isolated Dirac points. It is known that the third component of the Berry curvature vanishes everywhere except at the Dirac points, where it diverges to infinity. So, people do not seem to define the Chern number in gapless systems because the Berry curvature is not well-defined. However, they have no problem using the $A$ integral to calculate the Berry phase around the degeneracy. This phase has even been experimentally detected in graphene, Weyl, etc systems.
The numerical scheme requires shifting the discretized parameter space by a small factor to keep the Berry curvature from blowing up numerically. Due to this, when I plot the Berry curvature over parameter/k-space, I see the Berry curvature peaked sharply near the Dirac points (instead of diverging to infinity). However, my issue is that I don't understand why this finite peak correctly gives the expected Berry phase. Summing $F$ in the vicinity of this peak correctly gives the $\pi$ Berry phase expected around the degeneracy. Yet, this peak seems to be a consequence of the numerical method (i.e. shifting by a small factor).
Why does this ill-defined Berry curvature give the correct Berry phase around a degeneracy, seemingly as a consequence of the numerical scheme?