anchor.diff:  
       end function
   
  +    pure elemental function distance_haversine(p1_lon, p1_lat, p2_lon, p2_lat, radius) result(dist)
  +        !! Haversine distance formula giving the distance (meters) between points p1, p2
  +        !! on a sphere with the specified radius (by default = radius_earth in meters)
  +        real(dp), intent(in) :: p1_lon, p1_lat
  +            !! Lon/lat coordinates of input point
  +        real(dp), intent(in) :: p2_lon, p2_lat
  +            !! Lon/lat coordinates 
  +        real(dp), optional, intent(in) :: radius
  +            !! Sphere radius. The output distance will have the same units as radius. By default use radius_earth
  +        real(dp) :: dist
  +
  +        ! Try to avoid round-off error in single precision case by using double-precision here
  +        real(force_double) :: a, dlat, dlon, r_
  +
  +        r_ = radius_earth * 1.0_force_double
  +        if(present(radius)) r_ = radius * 1.0_force_double
  +
  +        dlat = (p2_lat - p1_lat) * 1.0_force_double * DEG2RAD
  +        dlon = (p2_lon - p1_lon) * 1.0_force_double * DEG2RAD
  +        a = sin(dlat*0.5_force_double)**2 + cos(p1_lat*1.0_force_double*DEG2RAD)*&
  +            cos(p2_lat*1.0_force_double*DEG2RAD)*sin(dlon*0.5_force_double)**2
  +        a = min(a, 1.0_force_double)
  +        dist = 2.0_dp * atan2(sqrt(a), sqrt(1.0_dp - a)) * r_
  +
  +    end function
  +
       subroutine test_spherical_mod()
           !! Unit tests
           real(dp):: lon1, lon2, lat1, lat2, area
  -        real(dp):: ans
  +        real(dp):: ans, val
           ! the test cases assume this radius
           real(dp), parameter:: test_radius = 6378137.0_dp
           real(dp), parameter:: radius_adjustment_factor = (test_radius/radius_earth)**2
  
positive.code:             print*, 'FAIL'
          end if
  
  
          !! Test 2
          ans = 3379082.32316_dp
          lon1 = 108.0_dp
  
negative.code:             print*, 'FAIL', area * radius_adjustment_factor, ans
          end if
  
          !! Test 1 of distance_haversine
          lon1 = 0.0_dp
          lon2 = 1.0_dp
          lat1 = 0.0_dp
          lat2 = 0.0_dp
          ans = 111319.490793_dp ! In R, " distHaversine(c(0, 0), c(1, 0)) "
          val = distance_haversine(lon1, lat1, lon2, lat2, test_radius)
          if(abs(val - ans) < 1.0e-5_dp) then
              print*, 'PASS'
          else
              print*, 'FAIL distance_haversine 1', ans, val
          end if
  
          !! Test 2 of distance_haversine
          lon1 = 0.0_dp
          lon2 = 1.0_dp
          lat1 = 0.0_dp
          lat2 = 1.0_dp
          ans = 157425.537108_dp ! In R, " distHaversine(c(0, 0), c(1, 1)) "
          val = distance_haversine(lon1, lat1, lon2, lat2, test_radius)
          if(abs(val - ans) < 1.0e-5_dp) then
              print*, 'PASS'
          else
              print*, 'FAIL distance_haversine 2', ans, val
          end if
  
          !! Test 3 of distance_haversine
          lon1 = -179.0_dp
          lon2 = 12.0_dp
          lat1 = 56.0_dp
          lat2 = -72.0_dp
          ans = 18184349.8625_dp  ! In R, " distHaversine(c(-179, 56), c(12, -72)) "
          val = distance_haversine(lon1, lat1, lon2, lat2, test_radius)
          if(abs(val - ans) < 1.0e-3_dp) then
              print*, 'PASS'
          else
              print*, 'FAIL distance_haversine 3', ans, val
          end if
  
          !! Test 4 of distance_haversine -- using longitudes well outside [-180, 360]
          lon1 = -179.0_dp - 360.0_dp * 3
          lon2 = 12.0_dp + 360.0_dp * 2
          lat1 = 56.0_dp
          lat2 = -72.0_dp
          ans = 18184349.8625_dp  ! In R, " distHaversine(c(-179, 56), c(12, -72)) "
          val = distance_haversine(lon1, lat1, lon2, lat2, test_radius)
          if(abs(val - ans) < 1.0e-3_dp) then
              print*, 'PASS'
          else
              print*, 'FAIL distance_haversine 4', ans, val
          end if
  
      end subroutine
  
  end module
  