Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 5 commits into from
Mar 25, 2025
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
11e178a
Add a new geographic distance calculation routine in w3servmd and use…
dahonegger Jan 2, 2025
b7eb7b3
Merge branch 'NOAA-EMC:develop' into failed_nearest_point_fix
thesser1 Feb 24, 2025
4e1dbfc
Merge branch 'NOAA-EMC:develop' into failed_nearest_point_fix
thesser1 Mar 12, 2025
0cd3fd6
Corrected variable description
dahonegger Mar 24, 2025
5ea8344
Spacing aesthetics
dahonegger Mar 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 95 additions & 0 deletions model/src/w3servmd.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ MODULE W3SERVMD
!/ 01-Mar-2016 : Added W3THRTN and W3XYRTN for post ( version 6.02 )
!/ processing rotated grid data
!/ 15-Jan-2021 : Added UV_TO_MAG_DIR routine ( version 7.12 )
!/ 02-Jan-2025 : Added DIST_HAVERSINE routine ( version 7.xx )
!/
!/ Copyright 2009-2012 National Weather Service (NWS),
!/ National Oceanic and Atmospheric Administration. All rights
Expand Down Expand Up @@ -511,6 +512,100 @@ REAL FUNCTION EJ5P ( F, ALFA, FP, YLN, SIGA, SIGB )
!/ End of NEXTLN ----------------------------------------------------- /
!/
END FUNCTION EJ5P
!/ ------------------------------------------------------------------- /

REAL FUNCTION DIST_HAVERSINE ( LON1, LAT1, LON2, LAT2 )
!/
!/ +-----------------------------------+
!/ | WAVEWATCH III USACE/ERDC |
!/ | D. A. Honegger |
!/ | FORTRAN 90 |
!/ | Last update : 02-Jan-2025 |
!/ +-----------------------------------+
!/
!/ 02-Jan-2025 : Creation ( version 7.xx )
!/
! 1. Purpose :
!
! Computes the haversine distance between two points on a sphere
!
! 2. Method
!
! Haversine Formula: R.W. Sinnott, "Virtues of the Haversine",
! Sky and Telescope, vol. 68, no. 2, 1984, p. 159
!
! 3. Parameters :
!
! Parameter list
!
! ----------------------------------------------------------------
! LON1 Real I Longitude of 1st point
! LAT1 Real I Latitude of 1st point
! LON2 Real I Longitude of 2nd point
! LAT2 Real I Latitude of 2nd point
! ----------------------------------------------------------------
!
! 4. Subroutines used :
!
! None.
!
! 5. Called by :
!
! WW3_BOUNC
!
! 6. Error messages :
!
! 7. Remarks :
!
! This function uses the haversine formula, which is robust to
! rounding errors when calculating short distances
! (less than 1 km).
!
! 8. Structure :
!
! See source code.
!
! 9. Switches :
!
! None.
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
USE CONSTANTS
! DERA: Degrees to Radians (PI/180)
! RADE: Radians to Degrees (180/PI)
!/
!/ ------------------------------------------------------------------- /
!/ 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

! Compute the haversine of the central angle
A = SIN(DLAT / 2.0)**2 + COS(LAT1 * DERA) * COS(LAT2 * DERA) * SIN(DLON / 2.0)**2

! Compute the angular distance (c), ensuring no precision issues
C = 2.0 * ATAN2(SQRT(A), SQRT(MAX(0.0, 1.0 - A)))

! Compute the spherical distance
DIST_HAVERSINE = RADE * C

RETURN
END FUNCTION DIST_HAVERSINE
!/ ------------------------------------------------------------------- /

!/ ------------------------------------------------------------------- /
REAL FUNCTION DIST_SPHERE ( lo1,la1,lo2,la2 )
!/
Expand Down
18 changes: 15 additions & 3 deletions model/src/ww3_bounc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,9 @@ PROGRAM W3BOUNC
!/ | WAVEWATCH III NOAA/NCEP |
!/ | F. Ardhuin |
!/ | M. Accensi |
!/ | D. A. Honegger |
!/ | FORTRAN 90 |
!/ | Last update : 21-Jul-2020 |
!/ | Last update : 02-Jan-2025 |
!/ +-----------------------------------+
!/
!/ 24-May-2013 : Adaptation from ww3_bound.ftn ( version 4.08 )
Expand All @@ -38,6 +39,8 @@ PROGRAM W3BOUNC
!/ 15-May-2018 : Add namelist feature ( version 6.05 )
!/ 04-May-2020 : Update spectral conversion ( version 7.11 )
!/ 21-Jul-2020 : Support rotated pole grid ( version 7.11 )
!/ 02-Jan-2025 : Change geographic distance method ( version 7.xx )
!/ 02-Jan-2025 : Add verbose=2 display output ( version 7.xx )
!/
!/
!/ Copyright 2012-2013 National Weather Service (NWS),
Expand Down Expand Up @@ -138,7 +141,7 @@ PROGRAM W3BOUNC
USE W3IOBCMD, ONLY: VERBPTBC, IDSTRBC
USE W3IOGRMD, ONLY: W3IOGR
USE W3TIMEMD
USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE, DIST_SPHERE
USE W3SERVMD, ONLY: ITRACE, NEXTLN, EXTCDE, DIST_HAVERSINE
#ifdef W3_RTD
USE W3SERVMD, ONLY: W3EQTOLL
#endif
Expand Down Expand Up @@ -534,6 +537,9 @@ PROGRAM W3BOUNC
CALL CHECK_ERR(IRET)
END IF

! Display the location of the ingested NetCDF file
IF (VERBOSE.GE.2) WRITE(NDSO,*) 'FILEID:',IP,'LON:',LONS(IP),'LAT:',LATS(IP)

! freq and dir variables
IRET=NF90_INQ_VARID(NCID(IP),"frequency",VARID(4))
CALL CHECK_ERR(IRET)
Expand Down Expand Up @@ -649,7 +655,7 @@ PROGRAM W3BOUNC
DO IP=1,NBO2
! Searches for the nearest 2 points where spectra are available
IF (FLAGLL) THEN
DIST=DIST_SPHERE ( LONS(IP),LATS(IP),XBPO(IP1),YBPO(IP1) )
DIST=DIST_HAVERSINE( LONS(IP),LATS(IP),XBPO(IP1),YBPO(IP1) )
ELSE
DIST=SQRT((LONS(IP)-XBPO(IP1))**2+(LATS(IP)-YBPO(IP1))**2)
END IF
Expand All @@ -671,6 +677,12 @@ PROGRAM W3BOUNC
END IF
END IF
END IF
! Display iteration to find distance minima
IF (VERBOSE.GE.2) WRITE(NDSO,*) &
'BOUNDID:',IP1,'FILEID:',IP, &
'DX:',LONS(IP)-XBPO(IP1),'DY:',LATS(IP)-YBPO(IP1), &
'DIST:',DIST,'DMIN:',DMIN,'DMIN2:',DMIN2

END DO ! IP1=1,NBO2
IF (VERBOSE.GE.1) WRITE(NDSO,*) 'DIST:',DMIN,DMIN2,IP1,IPBPO(IP1,1),IPBPO(IP1,2), &
LONS(IPBPO(IP1,1)),LONS(IPBPO(IP1,2)),XBPO(IP1), &
Expand Down