great_circle_dist Function

 real(ESMF_KIND_R8) function great_circle_dist( q1, q2, radius )
      real(ESMF_KIND_R8), intent(IN)           :: q1(2), q2(2)
      real(ESMF_KIND_R8), intent(IN), optional :: radius
 
      real(ESMF_KIND_R8) :: p1(2), p2(2)
      real(ESMF_KIND_R8) :: beta
      integer n

      do n=1,2
         p1(n) = q1(n)
         p2(n) = q2(n)
      enddo

      beta = asin( sqrt( sin((p1(2)-p2(2))/2.)**2 + cos(p1(2))*cos(p2(2))*   &
                         sin((p1(1)-p2(1))/2.)**2 ) ) * 2.

      if ( present(radius) ) then
           great_circle_dist = radius * beta
      else
           great_circle_dist = beta   ! Returns the angle
      endif

  end function great_circle_dist