Calculate ww3_bounc geographic distances with haversines and utilize verbose=2 output

[dahonegger]

Pull Request Summary

Bug fix for #1344

Description

When assigning NetCDF spectra to input boundary nodes, 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.

This fix replaces usage of DIST_SPHERE in ww3_bounc with a new function, DIST_HAVERSINE, that is more robust when distances are on the order of 1 km or less. Both DIST_SPHERE and DIST_HAVERSINE are defined in w3servmd.

This fix additionally includes the option for more detailed diagnostic output from ww3_bounc that is triggered with the VERBOSE=2 namelist value in ww3_bounc.nml.

There are expected changes to regression tests tp2.8 and tp2.17 owing to modifications to the weight assignments of the linearly interpolated boundary spectra.

Issue(s) addressed

• Fixes #1344

Commit Message

Calculate ww3_bounc geographic distances with haversines and utilize verbose=2 output

Check list

• Branch is up to date with the authoritative repository (NOAA-EMC) develop branch.
• Checked the checklist for a developer submitting to develop.
• If a version number update is required, checked the updating version number checklist.
• If a new feature was added, a regression test for testing the new feature is added.

Testing

mww3_test_03/./work_PR3_UNO_MPI_d2_c                     (4 files differ)
mww3_test_03/./work_PR3_UNO_MPI_e_c                     (1 files differ)
mww3_test_03/./work_PR2_UQ_MPI_e                     (1 files differ)
mww3_test_03/./work_PR3_UNO_MPI_d2                     (3 files differ)
mww3_test_03/./work_PR2_UNO_MPI_d2                     (6 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2                     (9 files differ)
mww3_test_03/./work_PR2_UQ_MPI_d2                     (6 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2_c                     (8 files differ)
mww3_test_03/./work_PR1_MPI_e                     (1 files differ)
mww3_test_03/./work_PR1_MPI_d2                     (5 files differ)
mww3_test_03/./work_PR3_UQ_MPI_e                     (1 files differ)
mww3_test_09/./work_MPI_ASCII                     (0 files differ)
ww3_tp2.17/./work_ma                     (1 files differ)
ww3_tp2.17/./work_mc1                     (1 files differ)
ww3_tp2.17/./work_mb                     (1 files differ)
ww3_tp2.17/./work_mc                     (1 files differ)
ww3_tp2.17/./work_ma1                     (1 files differ)
ww3_tp2.17/./work_a                     (1 files differ)
ww3_tp2.17/./work_c                     (1 files differ)
ww3_tp2.17/./work_b                     (1 files differ)
ww3_tp2.6/./work_ST4_ASCII                     (0 files differ)
ww3_tp2.8/./work_PR3_UQ                     (4 files differ)

• How were these changes tested? Via matrix
• Are the changes covered by regression tests? (If not, why? Do new tests need to be added?) Tests that use ww3_bounc with FLAGLL=.true. cover the changes, but the newly utilized VERBOSE=2 value in ww3_bounc.nml is not tested.
• Have the matrix regression tests been run (if yes, please note HPC and compiler)? Yes, on an ERDC HPC with Intel
• Please indicate the expected changes in the regression test output, (Note the list of known non-identical tests.) Expected changes to the regression tests include: modifications to the weight assignments of the linearly interpolated boundary spectra.
• Please provide the summary output of matrix.comp (matrix.Diff.txt, matrixCompFull.txt and matrixCompSummary.txt):

matrixCompSummary.txt
matrixDiff.txt
matrixCompFull.txt

[MatthewMasarik-NOAA]

Thank you for submitting this fix @dahonegger. We will start looking at it and I'll reply back tomorrow with the CM assigned and the timeline.

[mickaelaccensi]

Thank you David for this improvement !

[MatthewMasarik-NOAA]

@dahonegger I spoke with @JessicaMeixner-NOAA offline and we will be assigning a CM on Monday when @sbanihash is back from leave. We will be back in touch then. Thank you again for this submission.

[sbanihash]

Thank you @dahonegger for your work on this fix. I would like to request @kestonsmith-noaa to take a look at this PR and review the changes. Either @JessicaMeixner-NOAA or I would also take a look as a second reviewer as well.

[kestonsmith-noaa]

On the second point, checking coordinate range, I would prefer not to check coordinate range within haversine, but rather at input. Longitude and latitude out of geographic range won't present numerical problems within haversine and there are cases it may be useful to have longitudes outside [-180,180], for instance in regional domains spanning the international date line. Best, Keston

…

On Mon, Mar 24, 2025 at 8:06 AM Saeideh Banihashemi < ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In model/src/w3servmd.F90 <#1392 (comment)>: > + ! (less than 1 km). + ! + ! 8. Structure : + ! + ! See source code. + ! + ! 9. Switches : + ! + ! None. + ! + ! 10. Source code : + ! + !/ ------------------------------------------------------------------- / + USE CONSTANTS + ! DERA: Degrees to Radians (PI/180) + ! RADE: Radius of the Earth If this is intended to be Radius of the Earth , the constants module does not define it that way. Here is what RADE is in constants.F90: REAL, PARAMETER :: RADE = 180. / PI !< RADE Conversion factor from radians to degrees ------------------------------ In model/src/w3servmd.F90 <#1392 (comment)>: > + !/ ------------------------------------------------------------------- / + !/ Parameter list + !/ + REAL, INTENT(IN) :: LON1, LAT1, LON2, LAT2 + !/ + !/ ------------------------------------------------------------------- / + !/ Local parameters + !/ + REAL :: DLON, DLAT, A, C + !/ + !/ ------------------------------------------------------------------- / + !/ + + ! Compute differences in latitude and longitude + DLAT = (LAT2 - LAT1) * DERA + DLON = (LON2 - LON1) * DERA Should we make sure LAT1, LAT2 are within [-90, 90] and LON1, LON2 are within [-180, 180] before computation? — Reply to this email directly, view it on GitHub <#1392 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AZUY35CXAQZWAFQQH5NY5ET2V7YLNAVCNFSM6AAAAABY7EYMSWVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDOMJQGEZDIOBYHA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[sbanihash]

On the second point, checking coordinate range, I would prefer not to check coordinate range within haversine, but rather at input. Longitude and latitude out of geographic range won't present numerical problems within haversine and there are cases it may be useful to have longitudes outside [-180,180], for instance in regional domains spanning the international date line. Best, Keston

@kestonsmith-noaa Thank you for your input. I was thinking more about the DLON calculations (as opposed to a [0 360] coordinate), in order to make sure the haversine of the central angle computation is correct. have to re-think to be honest. maybe it's not even important. Looking at the fact that the trigonometric functions involved are periodic, I think both coordinates should work.

[sbanihash]

Changes look good. Tests have been passed on Hera with expected changes in ww3_tp2.8 and tp_2.17.

---

********************* non-identical cases ****************************

---

mww3_test_03/./work_PR2_UNO_MPI_e (1 files differ)
mww3_test_03/./work_PR2_UQ_MPI_d2 (16 files differ)
mww3_test_03/./work_PR3_UQ_MPI_e_c (1 files differ)
mww3_test_03/./work_PR2_UNO_MPI_d2 (10 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2_c (16 files differ)
mww3_test_03/./work_PR2_UQ_MPI_e (1 files differ)
mww3_test_03/./work_PR3_UNO_MPI_e_c (1 files differ)
mww3_test_03/./work_PR1_MPI_d2 (11 files differ)
mww3_test_03/./work_PR3_UNO_MPI_d2_c (17 files differ)
mww3_test_03/./work_PR3_UNO_MPI_d2 (18 files differ)
mww3_test_03/./work_PR3_UQ_MPI_d2 (15 files differ)
mww3_test_09/./work_MPI_ASCII (0 files differ)
ww3_tp2.10/./work_MPI_OMPH (7 files differ)
ww3_tp2.16/./work_MPI_OMPH (4 files differ)
ww3_tp2.17/./work_b (1 files differ)
ww3_tp2.17/./work_mb (1 files differ)
ww3_tp2.17/./work_mc (1 files differ)
ww3_tp2.17/./work_a (1 files differ)
ww3_tp2.17/./work_ma1 (1 files differ)
ww3_tp2.17/./work_c (1 files differ)
ww3_tp2.17/./work_ma (1 files differ)
ww3_tp2.17/./work_mc1 (1 files differ)
ww3_tp2.6/./work_ST4_ASCII (0 files differ)
ww3_tp2.8/./work_PR3_UQ (8 files differ)
ww3_ufs1.3/./work_a (2 files differ)

matrixCompSummary.txt
matrixCompFull.txt
matrixDiff.txt

[sbanihash]

@dahonegger Thank you for your work on this fix. I have approved your PR, just waiting for @kestonsmith-noaa to also confirm that the changes look good and I will then merge the PR. Thank for your patience.

[kestonsmith-noaa]

The changes looks fine to me. My only other suggestion would be to include the units of the return variable in the description, given our own internal confusion on this point and change of units within the function. Something like:

!  1. Purpose :
!
!     Computes the haversine distance between two points on a sphere in decimal degrees

[sbanihash]

Thank you @kestonsmith-noaa . That would be actually helpful. @dahonegger would you mind making this small change? I can merge the PR afterwards. Thanks!

[aliabdolali]

Hi @sbanihash
David is on AL until next week. You can either wait or add it to the header and merge.
Thanks for helping with the review.

[sbanihash]

Hi @sbanihash David is on AL until next week. You can either wait or add it to the header and merge. Thanks for helping with the review.

Thanks @aliabdolali. Done!

[dahonegger]

Thanks all for the constructive feedback on this PR -- I wholeheartedly agree with the suggested modifications.

And, thank you @aliabdolali for your help making sure @sbanihash was not left waiting for a response on Tuesday!