anchor.diff:              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
  
positive.code: 
      end function
  
      subroutine test_spherical_mod()
          !! Unit tests
          real(dp):: lon1, lon2, lat1, lat2, area
          real(dp):: ans
          ! 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
  
negative.code:             print*, 'FAIL'
          end if
  
          !! Test 2
          ans = 3379082.32316_dp
          lon1 = 108.0_dp
  