anchor.diff:        enddo
         s = s + da
   !      write(*,*) s
  -      
  +
   !  Test for in or out
         if (360._rs - abs(s) <= tol) then
            sph_poly_inout = .true.
         else
            sph_poly_inout = .false.
         endif
  -      
  +
         return
      end function sph_poly_inout
   !===============================================================================
  
positive.code: ! Computes the endpoint given a starting point lon,lat, azimuth and angular distance
  
       implicit none
       
       real(rs),intent(in)  :: lon1_in,lat1_in,az_in,delta_in
       real(rs),intent(out) :: lon2,lat2
       real(rs)             :: lon1,lat1,az,delta
       logical,optional,intent(in) :: degrees
       logical          :: using_deg,deg
       
       lon1=lon1_in ; lat1=lat1_in ; az=az_in ; delta=delta_in
       
       using_deg = present(degrees)
       if (using_deg) then
          if (degrees) deg = degrees
       else
          deg = .false.
       endif
       
       if (using_deg .and. deg) then
         if (delta > 180.) then
            write(*,*)'spherical_geometry: step: Error: distance must be less than 180 degrees.'
  
negative.code:    subroutine sph2cart(phi_in,theta_in,r,x,y,z,degrees)
  !  Returns the cartesian coordinates from spherical ones
  !  Theta is colatitude, phi is longitude and r is radius
  
        implicit none
  
        real(rs),intent(in)  :: theta_in,phi_in,r
        real(rs),intent(out) :: x,y,z
        real(rs)             :: theta,phi
        logical,optional,intent(in)  :: degrees
  
        if (present(degrees)) then
           if (degrees) then
              theta = theta_in * pi / 1.8d2
  