	subroutine est_wood(head,d1,dw,io_out,ierr)
c
c     subroutine for estimation of wood-anderson amplitudes
c
c     input:  
c             head   - station/channel header
c             d1     - data array
c             dw     - spectral array
c             io_out - unit number for writing
c     output:  
c             ierr   - error return code
c
c
c     include files for header
c
      include '/usr/local/include/qlib2.inc'
      include 'station_hdr.inc'
c
      real*4 d1(*)
      real*4 dt
      real*4 amaxmm, pmax, xmag
      integer*4 npts, io_out, ierr
      integer*4 imax
      complex*8 dw(*)
      character*3 chan
      character*4 stat
      character*10 string,ch1,ch2,ch3,ch4
      logical lexist
c
c     initialize error return
c
      ierr = 0
c
c     set variables for Wood-Anderson calculation
c
      ind = 3
      ind2 = 0
      fl = 0.0
      fh = 0.0
      np = 0
c
      stat  = head.station_id
      chan  = head.channel_id
      npts  = head.num_pts
      dt  = 1.0/head.sample_rate
c
c     remove instrument and reconvolve with Wood Anderson
c
      call coinst2(npts,dt,ind,ind2,d1,dw,np,fl,fh,ierr)
      if (ierr .ne. 0) then
        return
      end if
c
c     amaxmm is the WA amplitude in millimeters
c     pmax is the period of the maximum
c
      call amaxper(npts,dt,d1,amaxmm,pmax,imax)
      xmag = alog10(amaxmm/pmax)
c
c     check file existence, other fun stuff
c
      inquire(file='wood.wave',exist=lexist)
      if (lexist) then
        open(unit=io_out,file='wood.wave',fileopt='eof')
      else
        open(unit=io_out,file='wood.wave')
        write(io_out,'(a2)') '# '
        write(io_out,'(a2)') '# '
      end if
c
c     this is ugly and awkward - but I cannot think of a clever
c       way to reproduce Seistool output
c
      call klen(stat,kk)
      write(string,'(f10.3)') amaxmm
      call kstart(string,k1)
      call klen(string,k2)
      ch1 = string(k1:k2)
      id1 = k2 - k1 + 1
      write(string,'(f10.3)') pmax
      call kstart(string,k1)
      call klen(string,k2)
      ch2 = string(k1:k2)
      id2 = k2 - k1 + 1
      write(string,'(i10)') imax - 1
      call kstart(string,k1)
      call klen(string,k2)
      ch3 = string(k1:k2)
      id3 = k2 - k1 + 1
      write(string,'(f10.3)') xmag
      call kstart(string,k1)
      call klen(string,k2)
      ch4 = string(k1:k2)
      id4 = k2 - k1 + 1
c
c     this should be identical to a Seistool WA amplitude file
c
      write(io_out,100) stat(1:kk),chan,' WAS ',
     . ch1(1:id1),' ',ch2(1:id2),' ',ch3(1:id3),' ',
     . ch4(1:id4)
100   format(1x,10a)
      close(io_out)
c
c     all done
c
      return
      end
