poly_area() gives wrong values for area when the side of polygon crosses a pole

[nikizadehgfdl]

Describe the bug
The exchange grid tool uses a function poly_area(), which returns the area of a polygon on globe with given vertices, to calculate the (exchange) grid cell areas. This function gives an unreasonably large values when the side of such polygons crosses a pole. This is because the line inegral formula which the function is based on is not valid when such crossing happens and is bigger than what it should be by pi*R**2 (which explains why the returned area is almost half the area of the Earth).

This pole crossing happens mostly for stretched ATM grids (non-GCA) giving rise to tiling errors for those grids and holes in the Antarctic land_mask.

To Reproduce
Get a version of frenctools with more verbosity for the exchange grid creation tool:

bash
git clone -b nnz_NONGCA_strechedGridIssues https://gitlab.gfdl.noaa.gov/fre/fre-nctools.git fre-nctools_gitlab
cd fre-nctools_gitlab
git checkout adde039  #adds calculated atm_area to land_mask files
export CONFIG_SITE=site-configs/gfdl-ws/config.site 
source site-configs/gfdl-ws/env.sh 
module load gcc/6.2.0
export PATH=.:$PATH
autoreconf -i
./configure 
make

Now create a coupled mosaic with a stretched grid for ATM such as c256r25tlat32.0_om4p25.
Explore the values of area_atm array along the longitude in the tile that contains the south pole (land_mask_tile3.nc) . See figure below.

Expected behavior
area_atm should not give huge values.

System Environment
fre-nctools on PAN

Additional context
See issue #44 for the inaccuracies in poly_area even for non-stretched grids.

Figure below shows the value of area_atm in land_mask_tile3.nc

pathname320='/work/Niki.Zadeh/MOM6-examples_myfork/ice_ocean_SIS2/OM4_025/preprocessing.20201028/mosaic_c256r25tlat32.0_om4p25_githubMaster_Verbose_20210121/'          
land_mask0 = nc.Dataset(pathname320+'land_mask_tile3.nc')
mask320= np.array(land_mask0.variables['mask'])
area_atm320= np.array(land_mask0.variables['area_atm'])
plt.plot(area_atm320[128,:],marker='o')

Same as above after a fix to poly_area formulation to take into account the pole crossing for a polygon edge.

There is still the discontinuity mentioned in issue #44