diff --git a/model/src/w3servmd.F90 b/model/src/w3servmd.F90 index 600e6e572..7973d1c80 100644 --- a/model/src/w3servmd.F90 +++ b/model/src/w3servmd.F90 @@ -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 @@ -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: Radius of the Earth + !/ + !/ ------------------------------------------------------------------- / + !/ 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 ) !/ diff --git a/model/src/ww3_bounc.F90 b/model/src/ww3_bounc.F90 index 549f027d5..04956c7a3 100644 --- a/model/src/ww3_bounc.F90 +++ b/model/src/ww3_bounc.F90 @@ -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 ) @@ -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), @@ -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 @@ -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) @@ -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 @@ -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), &