anchor.diff:  !  Returns azimuth from point 1 to point 2.
   !  From: http://www.movable-type.co.uk/scripts/latlong.html
   
  -	  implicit none
  -
  -	  real(rs) :: azimuth,lon1,lat1,lon2,lat2
  -	  real(rs) :: rlon1,rlat1,rlon2,rlat2,d,dlon,dlat,conversion
  -	  logical,optional :: degrees
  -	  
  -	  conversion = 1._rs
  -	  if (present(degrees)) then
  -	     if (degrees) conversion = pi/180._rs
  +     implicit none
  +
  +     real(rs) :: azimuth,lon1,lat1,lon2,lat2
  +     real(rs) :: rlon1,rlat1,rlon2,rlat2,conversion
  +     logical,optional :: degrees
  +     
  +     conversion = 1._rs
  +     if (present(degrees)) then
  +        if (degrees) conversion = pi/180._rs
        endif
  -	  	  
  -	  rlon1 = conversion*lon1  ;  rlon2 = conversion*lon2
  -	  rlat1 = conversion*lat1  ;  rlat2 = conversion*lat2
  -	  
  +          
  +     rlon1 = conversion*lon1  ;  rlon2 = conversion*lon2
  +     rlat1 = conversion*lat1  ;  rlat2 = conversion*lat2
  +     
        azimuth = atan2(sin(rlon2-rlon1)*cos(rlat2) , &
                  cos(rlat1)*sin(rlat2) - sin(rlat1)*cos(rlat2)*cos(rlon2-rlon1) )
   
  -	  if (azimuth < 0) then
  -		  azimuth = azimuth+2._rs*pi
  -	  endif
  -	  
  -	  azimuth = azimuth / conversion
  -	  
  -!	  write(*,*)'Azimuth',azimuth
  -	  
  -	  return
  -	  
  +     if (azimuth < 0) then
  +        azimuth = azimuth+2._rs*pi
  +     endif
  +     
  +     azimuth = azimuth / conversion
  +     
  +!     write(*,*)'Azimuth',azimuth
  +     
  +     return
  +     
      end function azimuth
   !==============================================================================
   
  
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(*,*)'Error: distance must be less than 180 degrees.'
  			 stop
  		 else if ( lon1 <-180 .or. lon1 > 180 ) then
  			 write(*,*)'Error: longitude must be in range -180 - 180.'
  			 stop
  		 else if ( lat1 <-90 .or. lat1 > 90 ) then
  			 write(*,*)'Error: latitude must be in range -90 - 90.'
  			 stop
  		 endif
  	  else
  		 if (delta > pi) then
              write(*,*)'Error: distance must be less than 2pi radians.'
              stop
  	     else if (lon1 < -pi .or. lon2 > pi) then
              write(*,*)'Error: longitude must be in range -2pi - 2pi.'
              stop
  	     else if (lat1 < -pi/2.d0 .or. lat2 > pi/2.d0) then
              write(*,*)'Error: latitude must be in range -pi - pi.'
  		    stop
  	     endif
  	  endif
  	  
  	  if (using_deg .and. deg) then
  !  Convert to radians
  		  lon1=lon1*pi/1.8D2 ; lat1=lat1*pi/1.8D2
  		  az=az*pi/1.8D2     ; delta=delta*pi/1.8D2
  	  endif
  	  
  !  Calculate point which is delta degrees/radians from lon1,lat1 along az
       lat2 = asin(sin(lat1)*cos(delta) + cos(lat1)*sin(delta)*cos(az))
       lon2 = lon1 + atan2(sin(az)*sin(delta)*cos(lat1),  &
                           cos(delta)-sin(lat1)*sin(lat2) )
          
  	  if (using_deg .and. deg) then
  !  Convert to degrees
  		  lat2=1.8D2*lat2/pi  ; lon2=1.8D2*lon2/pi
  		  if(lon2>1.8D2) lon2=lon2-3.6D2 ; if(lon2<-1.8D2) lon2=lon2+3.6D2
  	  end if
  	  
  	  return
     
     end subroutine step
  !==============================================================================
  
negative.code: !  delta returns the angular distance between two points on a sphere given 
  !  the lat and lon of each using the Haversine formula
  !
       implicit none
       real(rs) :: delta,lat1,lon1,lat2,lon2
       logical,optional :: degrees
            
       if (present(degrees)) then
          if (degrees) then
           lat1=lat1*pi/1.8D2 ; lon1=lon1*pi/1.8D2
           lat2=lat2*pi/1.8D2 ; lon2=lon2*pi/1.8D2
          endif
       endif
       
       delta=atan2( sqrt( (cos(lat2)*sin(lon2-lon1))**2 + (cos(lat1)*sin(lat2) - &
             sin(lat1)*cos(lat2)*cos(lon2-lon1))**2) , &
             sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1))
  
       if (present(degrees)) then
           if (degrees) delta = delta * 1.8D2/pi
        endif
       
       return
       
     end function delta
  !==============================================================================
  
  