subroutine find_south_pole(rtime, glat, glon)

  real, intent(in)  :: rtime
  real, intent(out) :: glat, glon

  real    :: dl, alat, alon, a, bmag, xmag, ymag, zmag, v
  real    :: old
  integer :: i
  logical :: found, inner_found

  glat = -70.0

  found = .false.

  dl = 16.0

  glon = -dl

  old = 0.0

  do while (.not.found)
  
     call apex(rtime, glat, glon + dl, 110.0, a, &
          alat, alon, bmag, xmag, ymag, zmag, v)

     if (alat < old) then
        glon = glon + dl
        old  = alat 
     else
        if (dl < 0.5) then
           found = .true.
        else
           glon = glon - dl
           old  = 0.0
           dl   = dl / 2.0
        endif
     endif

  enddo

  found = .false.

  dl = 1.0

  old = 0.0

  do while (.not.found)
  
     call apex(rtime, glat, glon - dl, 110.0, a, &
          alat, alon, bmag, xmag, ymag, zmag, v)

     if (alat < old) then
        glat = glat - dl
        old  = alat 
     else
        if (dl < 0.5) then
           found = .true.
        else
           glat = glat + dl
           old  = 0.0
           dl   = dl / 2.0
        endif
     endif

  enddo

  return

end subroutine find_south_pole

subroutine find_geographic_point(rtime, alat, alon, glat, glon, IsVerbose)

  implicit none

  real, intent(in)  :: rtime, alat, alon
  logical, intent(in) :: IsVerbose
  real, intent(out) :: glat, glon

  real    :: dl, alattest, alontest, a, bmag, xmag, ymag, zmag, v
  real    :: old, glat0, glon0, glat_save, glon_save, alonin, alatin
  integer :: i, j, ifrom, ito
  logical :: found, inner_found, IsInitialized=.false.

  real :: dla, dlo, dist, disttest

  real, dimension(0:360,-90:90,2) :: SavedGeoToMagCoords

  save SavedGeoToMagCoords

  if (.not.IsInitialized) then
     IsInitialized = .true.
     SavedGeoToMagCoords = -100.0
  endif

  alonin = mod(alon+360.0,360.0)
  alatin = alat
  found = .false.

  if (IsVerbose) write(*,*) "Entering find_geographic_point"

  if (alat < -90 .or. alat > 90) then
     found = .false.
  else
     if (SavedGeoToMagCoords(int(alonin),int(alatin),1) > -100) then
        glat = SavedGeoToMagCoords(int(alonin),int(alatin),1)
        glon = SavedGeoToMagCoords(int(alonin),int(alatin),2)
        return
     endif
  endif

  alonin = alon

  if (alonin > 180.0) alonin = alonin - 360.0

  glon = alonin
  glat = alat

  found = .false.

  dla =  5.0
  dlo = 30.0

  do while (.not.found)

     dist = 1000.0

     glat0 = glat
     glon0 = glon

     ifrom = -2
     ito = 2

     if (ifrom*dla+glat0 < -90.0) ifrom = (-90.0-glat0)/dla
     if (ito*dla+glat0 > 90.0) ito = (90.0-glat0)/dla

     if (IsVerbose) write(*,*) glat, glon

     do i=ifrom,ito
        do j=-3,3

           glat = int(glat0 + i*dla)
           glon = int(glon0 + j*dlo)
           glon = mod(glon,360.0)

           if (IsVerbose) write(*,*) glat, glon

           call apex(rtime, glat, glon, 110.0, a, &
                alattest, alontest, bmag, xmag, ymag, zmag, v)

           disttest = abs(alattest - alat) + abs(alontest - alonin)

           if (disttest < dist) then
              glat_save = glat
              glon_save = glon
              dist = disttest
           endif

        enddo
     enddo

     glat = glat_save
     glon = glon_save

     call apex(rtime, glat, glon, 110.0, a, &
          alattest, alontest, bmag, xmag, ymag, zmag, v)

     dist = sqrt((alattest - alat)**2 + (alontest - alonin)**2)

     if (dist < 2.5) found = .true.
     if (dlo < 1.0) found = .true.

     dla = dla/2.0
     dlo = dlo/2.0

  enddo

  if (dist < 2.5 .and. alatin > -90 .and. alatin < 90) then
     alonin = mod(alon+360.0,360.0)
     SavedGeoToMagCoords(int(alonin),int(alatin),1) = glat
     SavedGeoToMagCoords(int(alonin),int(alatin),2) = glon
  endif

  return

end subroutine find_geographic_point