anchor.diff:  
   subroutine update_auto_stations()
   
  -use o92util, only: displacement_file, &
  +use earth, only: radius_earth_km
  +use geom, only: distaz2lola
  +use trig, only: d2r
  +
  +use o92util, only: coord_type, &
  +                   displacement_file, &
                      disp_file_save, &
                      station_file, &
                      nstations, &
                      stations, &
  +                   centroid, &
                      auto_n, &
  -                   auto_depth
  +                   auto_depth, &
  +                   auto_az, &
  +                   auto_mode
   
   implicit none
   
   ! Local variables
  -integer :: i, j
  -double precision :: disp(2002,6), xmin, xmax, ymin, ymax
  +integer :: i, j, ierr
  +double precision :: disp(2002,6), disp_mag, xmin, xmax, ymin, ymax
  +double precision :: d, lon, lat
   
   
   if (station_file.ne.'o92_autosta_86_this_when_finished') then
  
positive.code: 
  subroutine gcmdln()
  
  use io, only: stdout, verbosity
  use o92util, only: ffm_file, &
                     fsp_file, &
                     mag_file, &
  
negative.code:     moment = moment + faults(i,7)*faults(i,8)*faults(i,9)
  enddo
  centroid = centroid/moment
  write(stdout,1001) centroid(1),centroid(2),centroid(3)/1.0d3
  1001 format('auto_stations: centroid(lon,lat,dep(km)): ',2(F10.3),F9.2)
  
  
  if (auto_mode.eq.'h') then
  
      ! Automatic horizontal grid
  
      ! Sample along transects extending 500 km north, south, east, and west of centroid
      if (coord_type.eq.'geographic') then
          ! West-east
          dx = -500.0d0
          do while (dx.le.500.0d0)
              call distaz2lola(centroid(1),centroid(2),dx/radius_earth_km,90.0d0,lon,lat, &
                               'radians','degrees',ierr)
              if (ierr.ne.0) then
                  call usage('auto_stations: error computing longitude and latitude')
              endif
              write(81,*) lon,lat,auto_depth
              dx = dx + 1.0d0
          enddo
  
          ! South-north
          dy = -500.0d0
          do while (dy.le.500.0d0)
              call distaz2lola(centroid(1),centroid(2),dy/radius_earth_km,0.0d0,lon,lat, &
                               'radians','degrees',ierr)
              write(81,*) lon,lat,auto_depth
              if (ierr.ne.0) then
                  call usage('auto_stations: error computing longitude and latitude')
              endif
              dy = dy + 1.0d0
          enddo
  
      elseif (coord_type.eq.'cartesian') then
          ! Left-right (west-east)
          dx = -500.0d0
          do while (dx.le.500.0d0)
              x = centroid(1) + dx
              write(81,*) x,centroid(2),auto_depth
              dx = dx + 1.0d0
          enddo
  
          ! Bottom-top (south-north)
          dy = -500.0d0
          do while (dy.le.500.0d0)
              y = centroid(2) + dy
              write(81,*) centroid(1),y,auto_depth
              dy = dy + 1.0d0
          enddo
  
      endif
  
  
  elseif (auto_mode.eq.'v') then
  
      ! Automatic vertical grid
  
      ! Sample along transects extending 500 km positive/negative azimuth, up/down of centroid
      if (coord_type.eq.'geographic') then
          ! Along azimuth
          dx = -500.0d0
          do while (dx.le.500.0d0)
              call distaz2lola(centroid(1),centroid(2),dx/radius_earth_km,auto_az,lon,lat, &
                               'radians','degrees',ierr)
              if (ierr.ne.0) then
                  call usage('auto_stations: error computing longitude and latitude')
              endif
              write(81,*) lon,lat,centroid(3)/1.0d3
              dx = dx + 1.0d0
          enddo
      elseif (coord_type.eq.'cartesian') then
          ! Along azimuth
          dx = -500.0d0
          do while (dx.le.500.0d0)
              x = centroid(1) + dx*sin(auto_az*r2d)
              y = centroid(2) + dx*cos(auto_az*r2d)
              write(81,*) x,y,centroid(3)/1.0d3
              dx = dx + 1.0d0
          enddo
      endif
  
      ! Vertical
      dy = -500.0d0
      do while (dy.le.500.0d0)
          y = centroid(3)/1.0d3 + dy
          dy = dy + 1.0d0
          write(81,*) centroid(1),centroid(2),y
      enddo
  
  else
      call usage('auto_stations: no auto_mode named '//trim(auto_mode)//' (usage:station)')
  endif
  
  
  close(81)
  
  