anchor.diff:  C----
         a = sin(dlat/2.0e0)*sin(dlat/2.0e0) +
        1   cos(lat1*d2r)*cos(lat2*d2r)*sin(dlon/2.0e0)*sin(dlon/2.0e0)
  -      if (a.eq.1.0) then
  +      if (a.ge.1.0) then
             dist = 0.0e0
         else
             dist = 2.0e0*atan2(sqrt(a),sqrt(1.0-a))
  
positive.code: C----
        a = dsin(dlat/2.0d0)*dsin(dlat/2.0d0) +
       1   dcos(lat1*d2r)*dcos(lat2*d2r)*dsin(dlon/2.0d0)*dsin(dlon/2.0d0)
        if (a.eq.1.0d0) then
            dist = 0.0d0
        else
            dist = 2.0d0*datan2(dsqrt(a),dsqrt(1.0-a))
  
negative.code: C (Double precision)
  C----
        IMPLICIT none
        REAL*8 pi,d2r
        PARAMETER (pi=4.0d0*datan(1.0d0),d2r=pi/1.8d2)
        REAL*8 lon1,lat1,lon2,lat2,colat1,colat2,dlon,dlat,a
        REAL*8 dist,az
  C----
  C Check if points are polar opposite
  C----
        if (dabs(lat1+lat2).le.1.0d-6.and.
       1                  dabs(mod(lon1-lon2,1.8d2)).le.1.0d-6) then
            dist = pi
            az = 0.0d0
            goto 11
  