pure function great_circle_distance(r,long1,lat1,long2,lat2) result(d)

    implicit none

    real(wp)            :: d        !! great circle distance from 1 to 2 [km]
    real(wp),intent(in) :: r        !! radius of the body [km]
    real(wp),intent(in) :: long1    !! longitude of first site [rad]
    real(wp),intent(in) :: lat1     !! latitude of the first site [rad]
    real(wp),intent(in) :: long2    !! longitude of the second site [rad]
    real(wp),intent(in) :: lat2     !! latitude of the second site [rad]

    real(wp) :: c1,s1,c2,s2,dlon,clon,slon

    !Compute aux variables:
    c1    = cos(lat1)
    s1    = sin(lat1)
    c2    = cos(lat2)
    s2    = sin(lat2)
    dlon  = long1-long2
    clon  = cos(dlon)
    slon  = sin(dlon)

    d = r*atan2( sqrt((c2*slon)**2+(c1*s2-s1*c2*clon)**2), (s1*s2+c1*c2*clon) )

    end function great_circle_distance