great_circle Function

private pure function great_circle(r, lat1, lon1, lat2, lon2) result(d)

Great-circle distance (m) between two geographic points (deg), via the haversine formula (numerically stable for short arcs).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: r
real(kind=wp), intent(in) :: lat1
real(kind=wp), intent(in) :: lon1
real(kind=wp), intent(in) :: lat2
real(kind=wp), intent(in) :: lon2

Return Value real(kind=wp)


Called by

proc~~great_circle~~CalledByGraph proc~great_circle great_circle proc~spherical_tri_area spherical_tri_area proc~spherical_tri_area->proc~great_circle proc~tripolar_supergrid_arrays tripolar_supergrid_arrays proc~tripolar_supergrid_arrays->proc~great_circle proc~spherical_quad_area spherical_quad_area proc~tripolar_supergrid_arrays->proc~spherical_quad_area proc~metrics_fill_tripolar_whole metrics_fill_tripolar_whole proc~metrics_fill_tripolar_whole->proc~tripolar_supergrid_arrays proc~spherical_quad_area->proc~spherical_tri_area proc~metrics_fill_tripolar metrics_fill_tripolar proc~metrics_fill_tripolar->proc~metrics_fill_tripolar_whole proc~configure_ocean_metrics configure_ocean_metrics proc~configure_ocean_metrics->proc~metrics_fill_tripolar proc~engine_setup engine_setup proc~engine_setup->proc~configure_ocean_metrics

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: a
real(kind=wp), private :: dlam
real(kind=wp), private :: dphi
real(kind=wp), private :: h
real(kind=wp), private :: p1
real(kind=wp), private :: p2

Source Code

   pure function great_circle(r, lat1, lon1, lat2, lon2) result(d)
      !! Great-circle distance (m) between two geographic points (deg),
      !! via the haversine formula (numerically stable for short arcs).
      real(wp), intent(in) :: r, lat1, lon1, lat2, lon2
      real(wp) :: d
      real(wp) :: p1, p2, dphi, dlam, a, h
      p1 = lat1*DEG2RAD
      p2 = lat2*DEG2RAD
      dphi = (lat2 - lat1)*DEG2RAD
      dlam = (lon2 - lon1)*DEG2RAD
      ! wrap dlam into [-pi, pi] so the antimeridian doesn't blow up
      dlam = modulo(dlam + 3.14159265358979323846_wp, 2.0_wp*3.14159265358979323846_wp) &
             - 3.14159265358979323846_wp
      h = sin(0.5_wp*dphi)**2 + cos(p1)*cos(p2)*sin(0.5_wp*dlam)**2
      a = 2.0_wp*atan2(sqrt(h), sqrt(max(0.0_wp, 1.0_wp - h)))
      d = r*a
   end function great_circle