anchor.diff:              print*, 'FAIL'
           end if
   
  -
           !! Test 2
           ans = 3379082.32316_dp
           lon1 = 108.0_dp
  
positive.code:             print*, 'FAIL', area * radius_adjustment_factor, ans
          end if
  
      end subroutine
  
  end module
  
negative.code: 
      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, 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
  