	subroutine est_ground(head,d1,dw,rdist,io_out,iere)
c
c     subroutine for estimation of ground motion
c
c     input:  
c             stat   - station name
c             chan   - channel name
c             npts   - number of points in timeseries
c             dt     - sample spacing in sec
c             d1     - data array
c             dw     - spectral array
c             rdist  - distance from hypocenter to station in km
c                      =0 if not specified
c             io_out - unit number for writing
c     output:  
c             iere   - error return code
c
      parameter (maxdur = 4)
      parameter (maxp  = 65536)
c
c     include files for header and time
c
      include '/usr/local/include/qlib2.inc'
      record /INT_TIME/ t1_time
      record /EXT_TIME/ t2_time
	include 'station_hdr.inc'
c
      real*4 d1(*), d2(maxp)
      real*4 dt, hstep, totint
      real*4 amax_vel, pmax_vel
      real*4 amax_acc, pmax_acc
      real*4 amax_dsp, pmax_dsp
      integer*4 npts, io_out, iere
      double precision dt1
      integer*4 imax_vel, imax_acc, imax_dsp
      complex*8 dw(*)
      character*3 chan
      character*4 stat
      logical lexist
c
      real*4 gd(maxdur), sd(maxdur)
      integer*4 id(maxdur), npd(maxdur,2)
c
      data gd /0.05, 0.10, 0.20, 0.50/
c
c     acceleration of gravity in m/s/s
c
      grav = 9.81
c
c     initialize error return
c
      iere = 0
c
c     set variables for ground motion parameters 
c
      ind = 21
      ind2 = 6
      fl = 0.0
      fh = 0.1
	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 to velocity and hipass with a phaseless Butterworth
c
      call coinst2(npts,dt,ind,ind2,d1,dw,np,fl,fh,iere)
      if (iere .ne. 0) then
        return
      end if
c
c     find the velocity maximum (m/s) and its period
c
      call amaxper(npts,dt,d1,amax_vel,pmax_vel,imax_vel)
c
c     differentiate velocity to acceleration
c
      hstep = dt
      do kk = 1, npts - 1
        d2(kk+1) = (d1(kk+1) - d1(kk))/hstep
      end do
      d2(1) = 0.0
c
c     find the acceleration maximum (m/s/s) and its period
c
      call amaxper(npts-1,dt,d2,amax_acc,pmax_acc,imax_acc)
c
c     since we are here, determine the duration of strong shaking
c
      do kk = 1, maxdur
        id(kk) = 0
        npd(kk,1) = 0
        npd(kk,2) = 0
        do ii = 1, npts-1
          a = abs(d2(ii)/grav)
          if (a .ge. gd(kk)) then
            id(kk) = id(kk) + 1
            if (id(kk) .eq. 1) then
              npd(kk,1) = ii
            end if
            npd(kk,2) = ii
          end if
        end do
        if (id(kk) .ne. 0) then
          kpts = npd(kk,2) - npd(kk,1) + 1
          sd(kk) = kpts*dt
        else
          sd(kk) = 0.0
        end if
      end do
c
c     integrate velocity to displacement
c
      hstep = 0.5*dt
      totint = 0.0
      do kk = 1, npts - 1
        totint = totint + hstep*(d1(kk) + d1(kk+1))
        d2(kk) = totint
      end do
c
c     find the displacement maximum (m) and its period
c
      call amaxper(npts-1,dt,d2,amax_dsp,pmax_dsp,imax_dsp)
c
c     convert the measurements to millimeters from meters
c
      amax_dsp = amax_dsp * 1000.
      amax_vel = amax_vel * 1000.
      amax_acc = amax_acc * 1000.
c
c     output file for PGA,PGV,PGD
c
      call klen(stat,kk)
      inquire(file='ground.wave',exist=lexist)
      if (lexist) then
        open(unit=io_out,file='ground.wave',fileopt='eof')
      else
        open(unit=io_out,file='ground.wave')
        write(io_out,'(a2)') '# '
        write(io_out,'(a)') '#Stat Com  Nt  Val  Azi Dip  Year Doy Hr Mn
     . Sc    Ms         Amax    Freq(hz)    Lat      Lon     Dist(km)'
      end if
      dt1 = imax_acc * dble(dt)
      call f_add_dtime(head.begtime, dt1, t1_time)
      call f_int_to_ext(t1_time,t2_time)
      write(io_out,200) stat(1:kk), chan, head.network_id,' PGA ',
     .  head.azi, head.dip, t2_time.year, t2_time.doy, 
     .  t2_time.hour, t2_time.minute, t2_time.second, 
     .  int(t2_time.usec/usecs_per_tick),
     .  amax_acc, 1.0/pmax_acc, head.slat, head.slon, rdist
c
      dt1 = imax_vel * dble(dt)
      call f_add_dtime(head.begtime, dt1, t1_time)
      call f_int_to_ext(t1_time,t2_time)    
      write(io_out,200) stat(1:kk), chan, head.network_id,' PGV ',
     .  head.azi, head.dip, t2_time.year, t2_time.doy,
     .  t2_time.hour, t2_time.minute, t2_time.second,
     .  int(t2_time.usec/usecs_per_tick),
     .  amax_vel, 1.0/pmax_vel, head.slat, head.slon, rdist
c
      dt1 = imax_dsp * dble(dt)
      call f_add_dtime(head.begtime, dt1, t1_time)
      call f_int_to_ext(t1_time,t2_time)    
      write(io_out,200) stat(1:kk), chan, head.network_id,' PGD ',
     .  head.azi, head.dip, t2_time.year, t2_time.doy,
     .  t2_time.hour, t2_time.minute, t2_time.second,
     .  int(t2_time.usec/usecs_per_tick),
     .  amax_dsp, 1.0/pmax_dsp, head.slat, head.slon, rdist
c
200   format(1x,a4,1x,a3,1x,a4,a5,
     .  f4.0,1x,f4.0,1x,i4,1x,i3,1x,
     .  i2,1x,i2,1x,i2,1x,i5,1x,
     .  2f12.5,f7.3,1x,f8.3,1x,f12.5)
      close(io_out)
c
c     output file for shaking
c
      inquire(file='shake.wave',exist=lexist)
      if (lexist) then
        open(unit=io_out,file='shake.wave',fileopt='eof')
      else
        open(unit=io_out,file='shake.wave')
        write(io_out,'(a2)') '# '
        write(io_out,'(a)') '#Stat Com  Nt  Val  Azi Dip  Year Doy Hr Mn
     . Sc    Ms      5%    10%    25%    50%    Lat      Lon     Dist(km
     .)'
      end if
      write(io_out,250) stat(1:kk), chan, head.network_id,' SHAK',
     .  head.azi, head.dip, t2_time.year, t2_time.doy,
     .  t2_time.hour, t2_time.minute, t2_time.second,
     .  int(t2_time.usec/usecs_per_tick),
     .  sd(1), sd(2), sd(3), sd(4), head.slat, head.slon, rdist
c
250   format(1x,a4,1x,a3,1x,a4,a5,
     .  f4.0,1x,f4.0,1x,i4,1x,i3,1x,
     .  i2,1x,i2,1x,i2,1x,i5,1x,
     .  4f7.3,f7.3,1x,f8.3,1x,f12.5)
      close(io_out)
c
c     all done
c
      return
      end
