Add a new geographic distance calculation routine in w3servmd and use…

… this method in ww3_bounc to more robustly calculate distances between geographic locations. This commit also adds VERBOSE=2 diagnostic output. VERBOSE=2 was previously unused.

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 )
     !/

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,*) 'FILE-ID, LON, LAT: ',IP,LONS(IP),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), &