	subroutine est_energy(head,d1,dw,dist,iem,ic,io_out,ierr)
c
c     subroutine for estimation of energy magnitude
c
c     input:
c             head   - station/channel header
c             d1     - data array
c             dw     - spectral array
c             dist   - epicentral distance in km
c             iem    - crude component counter
c             ic     - component index
c             io_out - unit number for writing
c     output:
c             ierr   - error return code
c
      parameter (maxp  = 65536)
c
c     include files for header and time
c
      include '/usr/local/include/qlib2.inc'
      include 'station_hdr.inc'
c
      real*4 d1(*)
      real*4 data(maxp,3), d2(maxp)
      real*4 dt, dist, rdist, r0, r1, q0, q1, atten, atten2
      integer*4 npts, iem, io_out, ierr
      integer*4 ic(*)
      complex*8 dw(*)
      character*3 chan
      character*4 stat
      logical lexist
c
      external atten, atten2
c
c     initialize error flag
c
      ierr = 0
c
c     conversion from m to cm and from km to cm
c
      xm2cm = 100.
      xkm2m = 1000.
      xkm2cm = xkm2m * xm2cm
      pi = 4.0*atan(1.0)
c
c     reference depth and related parameters
c
      href = 8.0
      rho0 = 2.5
      beta0 = 3.0 * xkm2cm
      cf = 2.0
      fudge = 4.0 * pi * rho0 * beta0 / (cf**2)
c
c     station parameters
c
      stat  = head.station_id
      chan  = head.channel_id
      npts  = head.num_pts
      dt  = 1.0/head.sample_rate
      cid  = chan(3:3)
c
c     set variables for energy estimate
c
      ind = 21
      ind2 = 2
      fl = 0.05
      fh = 5.00
      np = 0
c
c     remove instrument to velocity and bandpass filter
c
      call coinst2(npts,dt,ind,ind2,d1,dw,np,fl,fh,ierr)
      if (ierr .ne. 0) then
        return
      end if    
c
c     save this data and convert to cm/s
c
      do ii = 1, npts
        data(ii,iem) = d1(ii) * xm2cm
      end do
c
c     if we don't have three-components yet, return
c
      if (iem .ne. 3) then
        return
      else
        do ii = 1, 3
          if (ic(ii) .ne. 1) then
            ierr = -iem
            return
          end if
        end do
      end if
c
c     form the array of velocities squared
c
      do kk = 1, npts
        d2(kk) = data(kk,1)*data(kk,1)
     .         + data(kk,2)*data(kk,2)
     .         + data(kk,3)*data(kk,3)
      end do
c
c     now form the total integral
c
      hstep = 0.5*dt
      totint = 0.0
      do kk = 1, npts - 1
        totint = totint + hstep*(d2(kk) + d2(kk+1))
      end do
c
c     now estimate the attenuation
c
      rdist = sqrt(dist*dist + href*href)
      r0 = href * xkm2cm
      r1 = rdist * xkm2cm
c
c     Either Richter attenuation, as modified by Jennings & Kanamori 1987
c
c      q0 = atten(href)
c      q1 = atten(rdist)
c
c     OR revised attenuation from Kanamori 1993
c
      q0 = atten2(href)
      q1 = atten2(rdist)
c
c     ok - form energy estimate
c
      atten_fac = ((r0*q0)/(r1*q1))**2
c
      ener = fudge * r1 * r1 * atten_fac * totint
      emag = (alog10(ener) - 9.05)/1.96
c
c     output file for energy magnitude
c
      call klen(stat,kk)
      inquire(file='energy.wave',exist=lexist)
      if (lexist) then
        open(unit=io_out,file='energy.wave',fileopt='eof')
      else
        open(unit=io_out,file='energy.wave')
        write(io_out,'(a2)') '# '
        write(io_out,'(a2)') '# '
      end if
      write(io_out,200) stat(1:kk),chan(1:2), head.network_id,
     .  ' EMAG ',emag, rdist, ener, totint
200   format(1x,a4,1x,a2,1x,a4,1x,a6,1x,f5.2,1x,f7.2,1x,2e15.5)
      close(io_out)
c
      do kk = 1, npts
        d1(kk) = d2(kk)
      end do
c
c     all done
c
      return
      end
