Non-unique nearest point calculation in ww3_bounc

[dahonegger]

Describe the bug
When assigning NetCDF spectra to input boundary nodes (MAPSTA=2, assigned via INBND_NML in ww3_grid.nml), ww3_bounc uses an iterative method to identify the two nearest neighbor spectra in spec.list. Custom diagnostic output in ww3_bounc reveals that points even one kilometer apart can be calculated as having a separation distance of zero. Moreover, some input point - boundary point pairs result in a returned distance of NaN. As a result, assignment of spectra to boundary nodes is sensitive to the ordering of the boundary nodes and to the ordering of the provided spectra in a manner that is not explicit in the code.

To Reproduce

1. Run ww3_grid
2. Provide multiple spectra that are within one kilometer of a boundary input node via files itemized in spec.list.
3. Run ww3_bounc

Expected behavior
Only one input spectrum is identified as having a minimum distance from each boundary node when using a grid/mesh with reasonable resolution.

Screenshot of the issue

Additional context
The screenshot consists of custom diagnostic output in ww3_bounc. Columns (labeled) are:

1. Boundary node index (outer loop)
2. Input spectrum file index (inner loop)
3. Distance in longitude between the boundary node and the input spectrum file
4. Distance in latitude between the boundary node and the input spectrum file
5. Total distance as calculated by DIST_SPHERE
6. Minimum distance from the boundary node to all input spectrum files (so far)
7. 2nd shortest distance from the boundary node to all input spectrum files (so far)

Note the presence of multiple locations with zero and NaN values for distance despite meridional and zonal distances (0.001-0.01 deg) that can reasonably be found along the boundaries of coastal grids/meshes. The example provided is at a location near the USACE FRF in Duck, North Carolina, USA (-75.7513, 36.1821).

Proposed fix
Calculate geographic distances with an updated formula that employs haversines. Haversine-based great-arc distances are more robust both to small angles and to floating point errors that can produce arccos arguments that are larger than one (and hence imaginary results). This fix, along with the additional lines of diagnostic stdout, is employed in a new branch in erdc/WW3, The screenshot below demonstrates the expected behavior.

Screenshot of the proposed fix