subroutine dist (lat1, lng1, lat2, lng2, dst) 

c  computes distance between two points on SPHERE 
c  given their latitudes and longitudes 
c  - returns arc distance, in degrees 

      real lat1,lng1,lat2,lng2,dst 

      rad = acos(-1.)/180. 

      dst = sin(lat1*rad)*sin(lat2*rad)+cos(lat1*rad)* 
     &      cos(lat2*rad)*cos((lng1-lng2)*rad) 

      dst = min(dst, 1.000000) 
      dst = max(dst,-1.000000) 
      dst = acos(dst)/rad 
      return 
      end