      program condet
c
c program to run sta/lta detection through continuous data stream,
c all input is read from condet.par
c
c two detection algorithms are supported:
c       - squared sta/lta
c       - Carl Johnson's 
c

c
c Lars Ottemoller, BGS, July 2005
c

c
c changes:
c
c 20/03/2008 lot - use several stations, network detection (option -net)
c    08/2008 lot - change write for extract.batch not to break up lines
c 14/08/2008 lot - added cosine taper to master signal when using COR
c                  this may give problem at file boundaries
c  2/10/2008 lot - fixed bug with wrong net detection at end and high ndet num
c 29/01/2009 lot - accept change of channel id
c 08/04/2010 lot - removed taper, and changed passes used in filter
c 30/04/2010 lot - added changes by WCC, mostly to netdet
c 15/09/2010 lot - fixed use of max_ratio
c  7/12/2010 lot - changed gt to ge in lines 1126 and 1135, as suggested by wc
c 14/12/2011 jh  - archive reading, replace reference to cwav and replace by wav
c 28/03/2012 jh  - space missing in extract.batch for network det.
c 19/12/2012 wcc - added corrections to lta and sta calculations
c 04/01/2013 jh  - remove reference to bud and scp, fix so work with arc
c 14/05/2013 lo  - merged last changes by wcc and jh
c  9/05/2014 pv  - fixed bug, condet could not read cor master file in arc mode
c 13 01 2016 jh  - removed subroutine hpsort

      implicit none 

c
c read SEISAN include files
c
      include 'seidim.inc'   ! dimensions
      include 'seisan.inc'   ! general seisan parameters
      include 'libsei.inc'
      include 'waveform.inc' ! for waveform handling
      include 'rea.inc'      ! for waveform handling
      integer seiclen
      character*4 base_type  ! type of data base like cont
      character*5 base_out   ! type for wavetool

      character*80 outfile   ! name of output seisan  waveform file
      character*80 text  
      character*1040 theader ! seisan channel header
      character*80 fheader(max_trace) ! seisan waveform main header
      character*29 headtext
      character*5 net_code   ! code of network
      real duration_min      ! minute
      integer i,j,l,k,sc     ! counter
      integer j_lta     	 ! counter WCC
      integer tl             ! number of samples to taper
      real taper,pi          ! taper value
      integer maxstat,maxtrig
      parameter (maxstat=100)
      character*5 station(maxstat),cbase(maxstat)
      character*4 component(maxstat)
      integer ind,ind_save
      double precision msec_start,msec_stop,msec,trigger_msec,msec1
      double precision triggerlen	! WCC
      double precision ptrigger_msec
      character*14 stop_time,start_time
      character*18 time
      integer year1,month1,day1,hour1,min1,isec1,doy1
      integer year2,month2,day2,hour2,min2,isec2,doy2
      real sec1,sec2
      real*8 sta,lta,star,ltar
      real*8 newlta,newsta	! WCC added to remove drift
      integer nsta,nlta
      real sta_duration,lta_duration
      logical startup
      integer write1,write2,write3,read1,write4
      logical b_flag
      integer code
      real trigger_ratio,detrigger_ratio,trigcor
      real trigger_min_int,trigger_min_dur
      logical trigger
      logical freeze
      character*8 filter_proto,filter_type  
      integer npoles,passes
      real flow,fhigh,samp_int
      integer ndc
      real rdc
      logical waveout
      character*80      data(max_data)  ! event array
      real data_save(max_sample),first_sta,first_lta
      real first_stacarl,first_ltacarl
      real ltacarl_save(max_sample),ltacarl(max_sample)
      real master(max_sample),cordata(max_sample)
      integer ntriggers(maxstat)
      character*10 detalg
      real carlratio,carlquiet
      integer nstat
      real pre_event,extduration
      logical netflag
      real netwindow,netminrat,netmaxdelt
      integer netmindet
      integer narg
      character*80 arg(40)
      character*80 master_filename
      integer nmaster ! number of samples in master signal
      real max
      integer cor_nsamp
      real max_ratio    ! WCC added maximum sta/lta ratio (for information only)

c
c   get seisan defaults including names of continuous waveform data bases
c   and archive data base
c
      call get_seisan_def 
c
c   get arguments
c
      call get_arguments(narg,arg)
      i=1
      netflag=.false.
      do while(i.le.narg)
        if (arg(i).eq.'-net') then
          netflag=.true.
        elseif (arg(i).eq.'-help') then
          write(*,*) ' usage: condet [-net] '
          write(*,*) '     -net: network detection using '//
     &      'condet.out as input '
          stop
        endif
        i=i+1
      enddo

      call get_contdet_def(nstat,base_type,cbase,start_time,stop_time,
     &   duration_min,station,component,waveout,flow,fhigh,
     &   sta_duration,lta_duration,trigger_ratio,
     &   detrigger_ratio,trigger_min_dur,trigger_min_int,
     &   freeze,detalg,
     &   carlratio,carlquiet,pre_event,extduration,
     &   netwindow,netmaxdelt,netminrat,netmindet,
     &   master_filename)
c
      base_out='-cwav'   ! default is to read from a a seisan cont
c
c      if(base_type(1:3).eq.'scp'.or.base_type(1:3).eq.'bud') then
c         arc=.true.   
c         arc_type=0   ! bud
c         base_out='-bud '
c         if(base_type.eq.'scp') then
c            arc_type=1
c            base_out='-scp '
c         endif
c      endif
      if(base_type(1:3).eq.'arc') then
        arc=.true.
        base_out='-arc '
      endif

      cseed=.false.

c
c open files
c
      if (.not.netflag) then
        call sei open( unknown$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'condet.out',    ! This filename.
     &               write1,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
        call sei open( unknown$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'nordic.out',        ! This filename.
     &               write2,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
        call sei open( unknown$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'condet.trace',    ! This filename.
     &               write4,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
      else
        call sei open( old$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'condet.out',    ! This filename.
     &               read1,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
        call sei open( unknown$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'net.out',    ! This filename.
     &               write1,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
        call sei open( unknown$,      ! Open unknown file.
     &               ' ',               ! No prompt.
     &               'condet.out.sort', ! This filename.
     &               write2,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
        call sei open( unknown$,      ! Open unknown file.
     &               ' ',               ! No prompt.
     &               'net.out.all',     ! This filename.
     &               write4,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)
      endif

      call sei open( unknown$,          ! Open old file.
     &               ' ',               ! No prompt.
     &               'extract.batch',   ! This filename.
     &               write3,            ! On this unit.
     &               b_flag,            ! Existance?.
     &               code )             ! Returned consition (n/a)

c
c call routine in case of network detection
c WCC
      if (netflag) then
        write(*,'(f8.1,1x,a)') extduration, ' extract duration'
                write(*,'(f8.1,1x,a)') pre_event,   ' pre event time'
        write(*,'(f8.1,1x,a)') netwindow,   ' net window sec'
        write(*,'(f8.1,1x,a)') netmaxdelt,  ' net max delta sec'
        write(*,'(f8.1,1x,a)') netminrat,   ' net min ratio'
        write(*,'(i8,1x,a)') netmindet,   ' net min detections'
        write(*,*) '--------------------------------------------'
        call netdet(read1,write1,write2,write3,write4,
     &     pre_event,extduration,
     &     netwindow,netmaxdelt,netminrat,netmindet,base_out)
        call sei close( close$, write1, code )
        call sei close( close$, write2, code )
        call sei close( close$, write3, code )
        call sei close( close$, write4, code )
        call sei close( close$, read1, code )
        stop
      endif


c
c init variables 
c
      npoles=4

c      if (cbase(1).eq.' ') then
      if (cbase(1).eq.' '.and..not.arc) then
        write(*,*) ' Name of continuous database '
        read(5,'(a)') cbase(1)
      endif

      if (start_time.eq.' ') then
        write(*,*) ' Start time (yyyymmdd...) '
        read(5,'(a)') start_time
      endif

      if (stop_time.eq.' ') then
        write(*,*) ' Stop time (yyyymmdd...) '
        read(5,'(a)') stop_time
      endif

      if (duration_min.eq.0.) then
        write(*,*) ' Duration in min '
        read(5,*) duration_min
      endif

      if (station(1).eq.' ') then
        write(*,*) ' Station '
        read(5,'(a)') station(1)
      endif

      if (component(1).eq.' ') then
        write(*,*) ' Component '
        read(5,'(a)') component(1)
      endif

c      write(*,*) ' start time: ',start_time
c      write(*,*) ' stop time:  ',stop_time
C WCC modified to add output of relevant parameters to stdout and condet.trace
       write(*,*) nstat,           ' stations'
       write(*,*) start_time,      ' start date'
       write(*,*) stop_time,       ' stop date'
       write(*,*) detalg,          ' detection algorithm'
       write(*,*) waveout,         ' waveout'
       write(*,*) duration_min,    ' interval to read (minutes)'
       write(*,*) flow,            ' filter low'
       write(*,*) fhigh,           ' filter high'
       write(*,*) trigger_min_int, ' min trigger interval (seconds)'
       write(*,*) extduration,     ' extract duration (seconds)'
       write(*,*) pre_event,       ' pre event time (seconds)'
       if (detalg(1:3).ne.'COR') then
          write(*,*) 'STA AND CARL ONLY PARAMETERS:'
          write(*,*) sta_duration,    '   sta duration'
          write(*,*) lta_duration,    '   lta duration'
          write(*,*) trigger_ratio,   '   trigger ratio'
          write(*,*) detrigger_ratio, '   detrigger ratio'
          write(*,*) trigger_min_dur, '   min trigger duration'
       else
          write(*,*) 'CORRELATION ONLY PARAMETERS:'
          write(*,*) master_filename, '   master waveform file'
          write(*,*) trigger_ratio,   '   trigger ratio (=100*corr min)'
       endif
       if (detalg(1:3).eq.'STA') then
          write(*,*) 'STA ONLY PARAMETERS:'
          write(*,*) freeze,          '   freeze lta'
       endif
       if (detalg(1:3).eq.'CAR') then
          write(*,*) 'CARL ONLY PARAMETERS:'
          write(*,*) carlratio,      '   carl ratio'
          write(*,*) carlquiet,      '   carl quiet'
       endif
c
c loop over stations
c
      do sc=1,nstat
c init
        ind_save=0
        nsta=0
        nlta=0
        sta=0.
        lta=0.
        star=0.
        ltar=0.
        trigger=.false.
        startup=.true.
        passes=1   ! clear recfil buffer when used first time
        filter_proto='BU'
        do k=1,max_sample
          data_save(k)=0.
          ltacarl_save(k)=0.
        enddo
c
        cont_base(1)=cbase(sc)
        ntriggers(sc)=0
        ptrigger_msec=0.
        trigger_msec=0.
c
c read master signal if detalg is correlation
c
        if (detalg(1:3).eq.'COR') then
          if (seiclen(master_filename).le.0) then
            write(*,*) ' MASTER WAVEFORM must be given in condet.par '
            stop
          endif
          call wav_init
          call get_full_wav_name(master_filename,wav_filename(1))
          call read_wav_header(1)
          ind=0
          do j=1,wav_nchan
            if (station(sc).eq.wav_stat(j).and.
     &          component(sc).eq.wav_comp(j).and.ind.eq.0) then
              ind=j
            endif
          enddo
          if (ind.eq.0) then
            write(*,*) ' channel not found ',station(sc),component(sc)
            stop
          else
            write(*,*) ' master file, channel selected index ',ind
          endif
          if(arc) then
            arc=.FALSE.                  ! turn arc off to read master file fron folder
            call wav_read_channel(ind)
            arc=.TRUE.                   ! turn arc back on
          else
            call wav_read_channel(ind)
          endif
          ndc=0
c put remove_dc back, lot 14/08/08
          call remove_dc(signal1,ndc,rdc,wav_nsamp(ind))
c this can be changed to compute dc and use it for taper, but not remove it
c as that is a problem when moving from one file to another
          do j=1,wav_nsamp(ind)
            write(15,*) j,signal1(j)
          enddo
c
c apply cosine taper 
c
          tl=int(wav_nsamp(ind)*.05)
          pi=acos(-1.)
          do j=1,tl
            taper=(-cos(pi/2.+(j-1)*(pi/2.)/float(tl-1)))
            signal1(j)=signal1(j)*taper
            signal1(wav_nsamp(ind)-j+1)=signal1(wav_nsamp(ind)-j+1)*
     &           taper
c            write(18,*) j,taper
          enddo

c          do j=1,wav_nsamp(ind)
c            write(16,*) signal1(j)
c          enddo
           
c
c filter data
c
          samp_int=1./wav_rate(ind)
c          if (flow.ne.0..or.fhigh.ne.0.) then
          if (flow.ne.0..and.fhigh.ne.0.) then
            filter_type='BP'
            call recfil(signal1,wav_nsamp(ind),signal1,
     &      filter_proto,0.,0.,npoles,filter_type,flow,fhigh,
     &      samp_int,passes)
          else
c
c apply highpass to remove dc
c
            filter_type='BP'   ! changed lo 16/09/2010
            flow=0.1
            fhigh=wav_rate(ind)/2.*.9
            call recfil(signal1,wav_nsamp(ind),signal1,
     &      filter_proto,0.,0.,npoles,filter_type,flow,fhigh,
     &      samp_int,passes)
          endif
c
c copy to master signal
c
          nmaster=wav_nsamp(ind)
          do j=1,wav_nsamp(ind)
            master(j)=signal1(j)
            write(17,*) master(j)
          enddo
        endif

c
c  signal that reading is from a continous data base
c
        cwav=.true.
        n_cont_base=1   ! for cont type
c
c   calculate end time and extended start time, needed in order to
c   make sure enough files read in
c
      cwav_start_time=start_time

      read(start_time,'(i4,5i2)') year1,month1,day1,hour1,min1,isec1
      sec1=float(isec1)
      call timsec(year1,month1,day1,hour1,min1,sec1,msec_start)

c 
c compute overall endtime
c
      read(stop_time,'(i4,5i2)') year2,month2,day2,hour2,min2,isec2
      sec2=float(isec2)
      call timsec(year2,month2,day2,hour2,min2,sec2,msec_stop)

c
c set initial values
c
      call cwav_time_limits(0)
      cwav_abs_start_time=msec_start
      cont_interval=duration_min*60.

      do while(cwav_abs_start_time.lt.msec_stop)

c
c set end of interval to be read
c
c        cwav_abs_end_time=cwav_abs_start_time+cont_interval
        call cwav_time_limits(1)
c
        if(.not.arc) then

c
c  read the header information for all files in cont bases in time
c  interval, assume info available in common block
c
            call cwav_read_bases
         else
c
c   read from archive
c
            call wav_read_arc_headers
c
c   stop if no data
c
            if(wav_nchan.eq.0) then
               write(6,*) ' No data, will stop'
               stop
            else
               write(6,*) ' Number of archive channels with data:',
     *         wav_nchan
            endif
          endif 
c
c find channel index
c
        ind=0

        do i=1,wav_nchan
          if (wav_stat(i).eq.station(sc).and.
     &        wav_comp(i).eq.component(sc)) 
     &    then
            ind=i
          endif
        enddo
        if (ind_save.eq.0.and.ind.ne.0) ind_save=ind

        if (ind_save.ne.ind) then
c          write(*,*) ' channel index changed, was/now ',
c     &        ind_save,ind
C WCC modified 2008/04/08 to allow channel index change
           write(*,*) station(sc),component(sc),
     &               ' channel index changed from ',
     &               ind, ' to ', ind_save
           if (ind.ne.0) ind_save=ind
        endif

c        cont_interval=duration_min*60.-
c     &       1./cwav_rate(ind,cwav_nseg(ind))

        cwav_abs_end_time=cwav_abs_start_time+cont_interval
 
        call sectim(cwav_abs_start_time,year1,doy1,month1,day1,hour1,
     &    min1,sec1)
        call sectim(cwav_abs_end_time,year2,doy2,month2,day2,hour2,
     &    min2,sec2)

        write(*,'(a,i4,"/",i2,"/",i2," ",i2,":",i2,":",f5.2)')
     &    ' start time: ',year1,month1,day1,hour1,min1,sec1
        write(write4,'(a,i4,"/",i2,"/",i2," ",i2,":",i2,":",f5.2)')
     &    ' start time: ',year1,month1,day1,hour1,min1,sec1
        write(*,'(a,i4,"/",i2,"/",i2," ",i2,":",i2,":",f5.2)')
     &    ' stop time:  ',year2,month2,day2,hour2,min2,sec2
        write(write4,'(a,i4,"/",i2,"/",i2," ",i2,":",i2,":",f5.2)')
     &    ' stop time:  ',year2,month2,day2,hour2,min2,sec2

        if (ind.ne.0) then
c          do j=1,cwav_nseg(ind)
c            write(*,*) ' segment:  ',cwav_nseg(j)
c            write(*,*) ' duration: ',cwav_duration(ind,j)
c            write(*,*) ' rate:     ',cwav_rate(ind,j)
c            write(*,*) ' samples:  ',cwav_nsamp(ind,j)
c          enddo

c
c  read the waveform data, one trace at a time
c
          call wav_read_channel(ind)    ! read trace
          if (wav_nsamp(ind).lt.0) wav_nsamp(ind)=0

          if (wav_nchan.eq.0) then
            startup=.true.
            goto 9999
          endif

c          write(*,*) ' total samples: ',wav_nsamp(ind)-1
          write(*,*) ' total samples: ',wav_nsamp(ind)
          write(write4,*) ' total samples: ',wav_nsamp(ind)
          write(*,*) ' channel index: ',ind
          write(write4,*) ' channel index: ',ind
C WCC modified 2009/04/08 to output component info
          write(*,*) ' ',station(sc),component(sc), 
     &                          ' channel index: ',ind
          write(write4,*) ' ',station(sc),component(sc), 
     &                          ' channel index: ',ind
c
c filter data
c

          samp_int=1./wav_rate(ind)
          if (flow.ne.0..and.fhigh.ne.0.) then
            filter_type='BP'
            call recfil(signal1,wav_nsamp(ind),signal1,
     &      filter_proto,0.,0.,npoles,filter_type,flow,fhigh,
     &      samp_int,passes)
          else
c
c apply highpass to remove dc
c
            filter_type='HP'
            flow=0.01
            call recfil(signal1,wav_nsamp(ind),signal1,
     &      filter_proto,0.,0.,npoles,filter_type,flow,fhigh,
     &      samp_int,passes)
          endif
          passes=-1   ! dont clear recfil buffer after used first time

c
c set waveform output variables
c
          if (waveout) then
            wav_out_nchan=2
            do j=1,wav_out_nchan
              wav_out_year(j)=wav_year(ind)
              wav_out_month(j)=wav_month(ind)
              wav_out_day(j)=wav_day(ind)
              wav_out_hour(j)=wav_hour(ind)
              wav_out_min(j)=wav_min(ind)
              wav_out_sec(j)=wav_sec(ind)
              wav_out_stat(j)=wav_stat(ind)
              if (j.eq.1) then
                wav_out_nsamp(j)=wav_nsamp(ind)
                wav_out_comp(j)=wav_comp(ind)
              else
                 wav_out_nsamp(j)=wav_nsamp(ind)
c                wav_out_nsamp(j)=wav_nsamp(ind)-1
                wav_out_comp(j)=detalg(1:4)
              endif
              wav_out_rate(j)=wav_rate(ind) 
              wav_out_cbyte(j)=wav_cbyte(ind)
            enddo
c
c   make and write seisan output header
c
            headtext=' '
c            net_code=detalg(1:3)
            net_code=wav_stat(ind)
            call sheads(wav_out_year,wav_out_month,
     *        wav_out_day,wav_out_hour,
     *        wav_out_min,wav_out_sec,wav_out_nchan,1,
     *        net_code,headtext,wav_out_stat,wav_out_comp,
     *        wav_out_nsamp,wav_out_rate,wav_out_cbyte,
     *        outfile,fheader,theader)
            outfile=outfile(1:seiclen(outfile))//'_'//detalg(1:3)
            write(6,*) ' Output file : '//outfile
            open(66,file=outfile,status='unknown',form='unformatted')

c
c   write headers, assume maximum of 30 channels, since only 12 headers are written out
c
            do i=1,12
              write(66) fheader(i)
            enddo
c
c write out filtered data
c
            call sheads(wav_out_year,wav_out_month,
     *        wav_out_day,wav_out_hour,
     *        wav_out_min,wav_out_sec,wav_out_nchan,1,
     *        wav_out_stat(1),headtext,wav_out_stat,wav_out_comp,
     *        wav_out_nsamp,wav_out_rate,wav_out_cbyte,
     *        text,fheader,theader)
            write(66) theader              
c find scale
            max=0.
            do i=1,wav_out_nsamp(1)
              if (abs(signal1(i)).gt.max) max=abs(signal1(i))
            enddo
            max=1000000./max
            write(66)(int(signal1(l)*max),l=1,wav_out_nsamp(1))      
c            write(66)(int(signal1(l))*100,l=1,wav_out_nsamp(1))      
          endif
c
c init rea
c
          call rea_hyp_clear(1)
          rea_id_line=' '
          if (waveout) then
            rea_nwav=1
            rea_wav(1)='                                        '//
     &               '                                       6'
            write(rea_wav(1)(2:79),'(a)') outfile(1:seiclen(outfile)) 
          endif

          hyp_year(1)=wav_year(1)
          hyp_month(1)=wav_month(1)
          hyp_day(1)=wav_day(1)
          hyp_hour(1)=wav_hour(1)
          hyp_min(1)=wav_min(1)
          hyp_sec(1)=wav_sec(1)
          hyp_dist_id(1)='L'
          rea_nphase=0
          rea_nhyp=1

          if (detalg(1:3).ne.'COR') then
            if (nsta.eq.0) then
                nsta=int(sta_duration*wav_rate(ind))
              write(*,*) ' nsta: ',nsta
            endif
            if (nlta.eq.0) then
              nlta=int(lta_duration*wav_rate(ind))
              write(*,*) ' nlta: ',nlta
            endif
            write(*,*) ' nsta/nlta: ',nsta,nlta
          endif
c
c compute cross correlation
c
          if (detalg(1:3).eq.'COR') then
            write(*,*) ' master nsamp ',nmaster,
     &        ' cont segment nsamp ',wav_nsamp(ind)
            call cor_time(master,nmaster,signal1,
c     &         cwav_nsamp(ind,cwav_nseg(ind)),
     &         wav_nsamp(ind),
     &         cordata,cor_nsamp)
c           do j=1,cwav_nsamp(ind,cwav_nseg(ind))
c              write(15,'(f5.2)') abs(cordata(j)*100.)
c           enddo
          endif

c
c compute sta and lta
c

          max_ratio=0.  ! WCC
          newlta=0.		! WCC
          newsta=0.		! WCC
          j_lta=0		! WCC

          do j=1,wav_nsamp(ind)   
          
C Important changes to the next two "if" constructs (WCC 7/2011)
C   1: Removed "+1" from indexing (made nsta and nlta one too short)
C   2: Changed (j.ge...) to (j.gt....) (avoid index 0 with the +1 removed)

c
c set previous data values for sta and lta
c
            if (detalg(1:3).ne.'COR') then
              if (j.gt.nsta) then
                first_sta=signal1(j-nsta)
              else
                first_sta=data_save(nlta-nsta+j)
              endif
              if (j.gt.nlta) then
                first_lta=signal1(j-nlta)
              else
                first_lta=data_save(j)
              endif
            endif
c
c previous values for carltrig
c
            if (detalg.eq.'CAR') then
              if (j.gt.nsta) then
                first_stacarl=ltacarl(j-nsta)
              else
                first_stacarl=ltacarl_save(nlta-nsta+j)
              endif
              if (j.gt.nlta) then
                first_ltacarl=ltacarl(j-nlta)
              else
                first_ltacarl=ltacarl_save(j)
              endif
            endif
c
c get sample time
c
            msec=cwav_abs_start_time+
     &         float(j-1)/wav_rate(ind)
c
c add new sample and delete first sample in sta and lta
c
            if (detalg(1:3).eq.'STA') then
              sta=sta+(signal1(j)**2 - first_sta**2)/float(nsta)
              ! NEW section calculates and synchronizes to sta at start of record
              if (j.lt.nsta) then							! WCC
              	newsta=newsta+signal1(j)**2/float(nsta)		! WCC
              elseif (j.eq.nsta) then						! WCC
              	sta=newsta+signal1(j)**2/float(nsta)		! WCC
              endif											! WCC
              ! End new section !!!!!!!!!!!!!!!!!!!!!!!!!!!!! WCC
              if (.not.trigger.or..not.freeze) then
                lta=lta+(signal1(j)**2 - first_lta**2)/float(nlta)
              	! NEW section calculates and synchronizes to lta at start of record
              	j_lta=j_lta+1								! WCC
              	if (j_lta.lt.nlta) then						! WCC
              		newlta=newlta+signal1(j)**2/float(nlta)	! WCC
              	elseif (j_lta.eq.nlta) then					! WCC
              		lta=newlta+signal1(j)**2/float(nlta) 	! WCC
                endif										! WCC
                ! End new section !!!!!!!!!!!!!!!!!!!!!!!!!!! WCC
              endif

              signal2(j)=sta/lta

            elseif (detalg(1:3).eq.'CAR') then
              sta=sta+signal1(j)/float(nsta)-
     &          first_sta/float(nsta)
              ltacarl(j)=ltacarl(j)+signal1(j)/float(nlta)-
     &            first_lta/float(nlta)
              star=star+abs(signal1(j)-ltacarl(j))/float(nsta)-
     &            abs(first_sta-first_stacarl)/float(nsta)
     
              ltar=ltar+abs(signal1(j)-ltacarl(j))/float(nlta)-
     &            abs(first_lta-first_ltacarl)/float(nlta)

              signal2(j)=
     &          star-carlratio*ltar
     &          -abs(sta-lta)
     &          -carlquiet

            elseif (detalg(1:3).eq.'COR') then
              signal2(j)=cordata(j)*100.
            endif

            if (signal2(j).gt.max_ratio) max_ratio=signal2(j)

            if (startup.and.j.ge.nlta) then
              startup=.false.
            endif
c
c save signal 
c
            if (detalg(1:3).ne.'COR') then
              if (j.eq.wav_nsamp(ind)) then
                do k=1,nlta
                  data_save(k)=signal1(wav_nsamp(ind)-nlta+k)		! WCC removed -1
                  ltacarl_save(k)=ltacarl(wav_nsamp(ind)-nlta+k)	! WCC removed -1
                enddo
              endif
            endif
c
c check for trigger based on sta/lta ratio
c
            if (.not.startup) then 
              if (signal2(j).gt.trigger_ratio.and.
     &        .not.trigger) then
                trigger_msec=msec
                trigger=.true.
                if (detalg.eq.'COR') trigcor=signal2(j)
              endif
              if (trigger) then
c
c check for detrigger, detection if trigger duration long enough
c
                if (signal2(j).le.detrigger_ratio) then
c                  write(*,*) ' detrigger ',signal2(j)

                  if (msec-trigger_msec.lt.trigger_min_dur) then
c                    write(*,*) ' trigger interval too short '
                    trigger=.false.
                  endif
                  if (trigger_msec-ptrigger_msec.lt.
     &              trigger_min_int) then
                      trigger=.false.
                  endif

                  if (detalg.eq.'COR'.and.trigcor.lt.signal2(j)) then
                    trigcor=signal2(j)
                    trigger_msec=msec
                  endif
                  if (trigger) then
c
c write out trigger
c
                    ntriggers(sc)=ntriggers(sc)+1

c
c add phase to nordic.out
c
                    rea_nphase=rea_nphase+1
                    if (rea_nphase.ge.max_phase) then
                      write(*,*) ' too many phases -> stop '
                      stop
                    endif
                    call rea_phase_clear(rea_nphase)
                    rea_stat(rea_nphase)=wav_stat(ind)
                    rea_comp(rea_nphase)=wav_comp(ind)
c                  rea_co(rea_nphase)(1:2)=wav_comp(ind)(1:1)//
c     &                 wav_comp(ind)(4:4)
                    rea_co(rea_nphase)(1:2)=wav_out_comp(2)(1:1)//
     &                 wav_out_comp(2)(4:4)
                    rea_phase(rea_nphase)=detalg(1:1)
                    rea_onset(rea_nphase)=' '
                    rea_weight_out(rea_nphase)='  '
                    rea_weight_in(rea_nphase)=' '
                    rea_polarity(rea_nphase)=' '
                    call sectim(trigger_msec,rea_year(rea_nphase),
     &     doy1,rea_month(rea_nphase),rea_day(rea_nphase),
     &     rea_hour(rea_nphase),rea_min(rea_nphase),
     &     rea_sec(rea_nphase))
  
                    call sectim(trigger_msec,year1,
     &                doy1,month1,day1,hour1,min1,sec1)
                    isec1=int(sec1)
                    time=' '
                    write(time,'(i4.4,5i2.2)') 
     &                year1,month1,day1,hour1,min1,isec1
                    if (detalg(1:3).ne.'COR') then
                      if (max_ratio.gt.9999.) max_ratio=9999.			! WCC
                      triggerlen=msec-trigger_msec						! WCC
                      if (triggerlen.gt.9999.9) triggerlen=9999.9		! WCC
                      write(write1,'(a5,1x,a4,1x,a14,1x,f6.1,1x,f6.1)') ! WCC
     &                  station(sc),component(sc),time,
     &                  triggerlen, max_ratio
         			  max_ratio=0.  ! WCC
                    else
                      write(write1,'(a5,1x,a4,1x,a14,1x,f6.1,1x,f6.1)')
     &                station(sc),component(sc),time,
     &                trigcor,0.  ! WCC
                    endif

c
c write out extract command
c
c                    msec1=trigger_msec-60.
                    msec1=trigger_msec-pre_event
                    call sectim(msec1,year1,
     &                doy1,month1,day1,hour1,min1,sec1)
                    isec1=int(sec1)
                    write(time,'(i4.4,4i2.2,f6.3)') 
c     &                year1,month1,day1,hour1,min1,isec1
     &                year1,month1,day1,hour1,min1,sec1
                     do k=1,18
                       if (time(k:k).eq.' ') time(k:k)='0'
                     enddo
c                    write(write3,*) 'wavetool -cwav -start '
                    write(write3,'(a,a,a,f6.1,a)')
     &                'wavetool ',base_out,' -start '
     &                // time(1:seiclen(time)) 
     &     // ' -duration ',extduration,' -wav_out_file SEISAN'
c
c save time as previous trigger
c
                    ptrigger_msec=trigger_msec
                    trigger=.false.
                  endif
                endif
              endif
            endif
          enddo
   
          if (waveout) then
c
c write out sta/lta
c
            call sheads(wav_out_year,wav_out_month,
     *        wav_out_day,wav_out_hour,
     *        wav_out_min,wav_out_sec,wav_out_nchan,2,
     *        wav_out_stat(2),headtext,wav_out_stat,wav_out_comp,
     *        wav_out_nsamp,wav_out_rate,wav_out_cbyte,
     *        text,fheader,theader)
            write(66) theader              
            write(66)(int(signal2(l)),l=1,wav_out_nsamp(2))      
          endif
c
c write out phases to Nordic file
c
          call rea_event_out(write2,.true.,data,code)
        endif
9999    continue
c
c set start of next interval to read
c
        cwav_abs_start_time=cwav_abs_end_time
      enddo
      enddo ! end of loop over stations
c
c close files
c
      call sei close( close$, write1, code )
      call sei close( close$, write2, code )
      call sei close( close$, write3, code )
      call sei close( close$, write4, code )
      do sc=1,nstat
        write(*,*) ' number of triggers ',station(sc),ntriggers(sc)
      enddo

      stop
      end


      subroutine get_contdet_def(nstat,base_type,
     *   cbase,startdate,stopdate,
     &   duration,station,component,waveout,flow,fhigh,
     &   sta_duration,lta_duration,trigger_ratio,
     &   detrigger_ratio,trigger_min_dur,trigger_min_int,
     &   freeze,detalg,
     &   carlratio,carlquiet,pre_event,extduration,
     &   netwindow,netmaxdelt,netminrat,netmindet,
     &   master_filename)

      implicit none
      include 'libsei.inc'
      include 'seidim.inc'

      character*80 def_file         ! file names
      integer in                    ! file unit
      integer code                  ! return code from opening def file
      integer i                     ! counter
      character*80 line             ! text line
      real var                      ! value of variables
      integer nstat                 ! station counter
      integer seiclen
      character*4 base_type       ! type of data base, archive type or seisan cont

      character*(*) cbase(*),startdate,stopdate,
     &  station(*),component(*),
     &  detalg
      logical waveout,freeze
      real duration,flow,fhigh,sta_duration,lta_duration
      real trigger_ratio,detrigger_ratio
      real trigger_min_int,trigger_min_dur
      real carlratio,carlquiet
      real pre_event,extduration
      real netwindow,netminrat,netmaxdelt
      integer netmindet
      character*80 master_filename
      real cormin
   
c
c check for file end open
c

      def_file = 'condet.par'

      call sei get file( open$+ignore$,   ! Open waveform file.
     &                   in,              ! On unit.
     &                   code,            ! Returned condition.
     &                   'DAT',           ! Alternative search directory.
     &                   def_file )       ! For this filename.

      
      if(code.ne.e_ok$) then
        write(*,*) ' definition file does not exist: condet.par'
        stop
      endif
c
c init values
c
      cbase(1)=' '
      base_type=' '
      startdate=' '
      stopdate=' '
      station(1)=' '
      component(1)=' '
      waveout=.false.
      duration=0.
      sta_duration=0.
      lta_duration=0.
      trigger_ratio=10.
      detrigger_ratio=10.
      trigger_min_int=0.
      trigger_min_dur=10.
      freeze=.false.
      detalg='STA'
      carlratio=4.
      carlquiet=2.
      nstat=0  ! station counter
      pre_event=60.
      extduration=180.
      netwindow=0.
      netmindet=1
      netminrat=0.
      netmaxdelt=0.
      cormin=1.
      master_filename=' '
      flow=0.
      fhigh=0.


100   continue

c
c read text line from file, check for code and set variables
c
      read(in,'(a80)',end=300,err=200) line

      if (line(1:10).eq.'START DATE') then
        read(line(41:54),'(a)') startdate
      elseif (line(1:9).eq.'BASE TYPE') then
        read(line(41:54),'(a)') base_type
      elseif (line(1:9).eq.'STOP DATE') then
        read(line(41:54),'(a)') stopdate
      elseif (line(1:7).eq.'STATION') then
        nstat=nstat+1
        read(line(41:45),'(a)') cbase(nstat)
        read(line(51:55),'(a)') station(nstat)
        read(line(57:60),'(a)') component(nstat)
      elseif (line(1:13).eq.'DET ALGORITHM') then
        read(line(41:50),'(a)') detalg
      elseif (line(1:15).eq.'MASTER WAVEFORM') then
        read(line(41:seiclen(line)),'(a)') master_filename
      elseif (line(1:7).eq.'WAVEOUT') then
        read(line(41:50),'(f10.1)') var
        if (var.eq.1.) waveout=.true.
      elseif (line(1:8).eq.'INTERVAL') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) duration=var
      elseif (line(1:15).eq.'CORRELATION MIN') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) cormin=var
      elseif (line(1:10).eq.'FILTER LOW') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) flow=var
      elseif (line(1:11).eq.'FILTER HIGH') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) fhigh=var
      elseif (line(1:10).eq.'STA LENGTH') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) sta_duration=var
      elseif (line(1:10).eq.'LTA LENGTH') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) lta_duration=var
      elseif (line(1:13).eq.'TRIGGER RATIO') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) trigger_ratio=var
      elseif (line(1:15).eq.'DETRIGGER RATIO') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) detrigger_ratio=var
      elseif (line(1:10).eq.'FREEZE LTA') then
        read(line(41:50),'(f10.1)') var
        if (var.eq.1.) freeze=.true.
      elseif (line(1:17).eq.'MIN TRIG DURATION') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) trigger_min_dur=var
      elseif (line(1:17).eq.'MIN TRIG INTERVAL') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) trigger_min_int=var
      elseif (line(1:10).eq.'CARL RATIO') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) carlratio=var
      elseif (line(1:10).eq.'CARL QUIET') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) carlquiet=var
      elseif (line(1:14).eq.'PRE EVENT TIME') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) pre_event=var
      elseif (line(1:14).eq.'NET WINDOW SEC') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) netwindow=var
      elseif (line(1:15).eq.'NET MAXDELT SEC') then             ! WCC
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) netmaxdelt=var
      elseif (line(1:13).eq.'NET MIN RATIO') then               ! WCC
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) netminrat=var
      elseif (line(1:11).eq.'NET MIN DET') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) netmindet=int(var)
      elseif (line(1:16).eq.'EXTRACT DURATION') then
        read(line(41:50),'(f10.1)') var
        if (var.gt.0.) extduration=var
      endif

      goto 100

200   continue
      write(*,*) ' Error in condet.par file'
      return

300   continue
      call sei close( close$, in, code )

      if (netmaxdelt.eq.0.) netmaxdelt=netwindow    ! WCC 3/2010

      if (detalg(1:3).eq.'CAR') then
        trigger_ratio=0.
        detrigger_ratio=0.
        sta_duration=1.0
        lta_duration=8.0
      elseif (detalg.eq.'COR') then
        trigger_min_dur=0.
        trigger_ratio=cormin*100.
        detrigger_ratio=0.
        sta_duration=0.
        lta_duration=0.
      endif

      return
      end


      subroutine netdet(read1,write1,write2,write3,write4,
     &     pre_event,extduration,netwindow,netmaxdelt,
     &     netminrat,netmindet,base_out)

      implicit none
      character*5 base_out
      integer read1,write1,write2,write3,write4	! WCC added 2&4
      real pre_event,extduration
      real netwindow, netmaxdelt, netminrat
      integer netmindet
      integer i,j,k,n,ind
      integer max
      parameter (max=500000)
      double precision msec(max),hmsec
      integer flag(max)	! WCC replaces used: 0=unused,1=used,2=rejected
      logical rpt
C      logical used(max)
      character*5 stat(max),hstat
      character*4 comp(max)	! WCC
      real coda(max),ratio(max)	! WCC
      integer year,month,day,hour,min,sec,doy
      real rsec
      integer ndet
	  integer ncorrdet	! WCC 3/2010
      character*14 time
      character*80 line

c
c read data    
c
      n=1
10    continue
      read(read1,'(a80)',end=50) line
c      write(*,*) line
      read(line,'(a5,1x,a4,1x,i4,5i2,1x,f6.1,1x,f6.1)') 
     &   stat(n),comp(n),year,month,day,hour,min,sec,coda(n),ratio(n)
      rsec=float(sec)
      call timsec(year,month,day,hour,min,rsec,msec(n))
      flag(n)=0
      n=n+1
      if (n.gt.max) then
      	write(*,*) 'More than ', max,
     &			' , only processing first ', max
      	write(write1,*) 'More than ', max,
     &			' detections, only processing first ', max
     	goto 50
      endif
      goto 10

50    continue
      n=n-1
      write(*,'(i8,1x,a)') n,' detections on all stations'
      write(write1,'(i8,1x,a)') n,' detections on all stations'

c WCC: Use heapsort, much faster
	call hpsort(n,msec,stat,comp,coda,ratio)
c
c now check if enough different stations in time window
c
      ndet=0
      ncorrdet=0	! WCC 3/2010
      do i=1,n
          hmsec=msec(i)
          call sectim(hmsec,year,
     &             doy,month,day,hour,min,rsec)
	      sec=int(rsec)
		  write (write2,'(a5,1x,a4,1x, i4,5i2.2,1x,f6.1,1x,f6.1)')
     &                 stat(i),comp(i),
     &                 year,month,day,hour,min,sec,coda(i),ratio(i)
          if (flag(i).eq.0 .and. ratio(i).ge.netminrat) then
              ndet=1
              j=i+1
              if (i.lt.n) then
                  do while(msec(j).le.msec(i)+netwindow.and.
     &                   msec(j).le.msec(j-1)+netmaxdelt.and.	! WCC 3/2010
     &                   j.le.n) ! lot 02/10/2008
                      rpt=.false.
C WCC Check that its trigger ratio is big enough
                if (ratio(j).ge.netminrat) then
C WCC verify that it's not an already-picked station
						  do k=i,j-1
							if (stat(j).eq.stat(k)) then
								rpt=.true.
							endif
						  enddo
						  if (.not.rpt) then
							  ndet=ndet+1
							  flag(j)=1
							  flag(i)=1
						  else
							  flag(j)=2		! Rejected
						  endif
				      else
				      	flag(j)=2	! Rejected
				      endif
                      j=j+1
                  enddo
              endif
			  if (ndet.ge.netmindet) then
				  hmsec=msec(i)-pre_event
				  call sectim(hmsec,year,
     &             doy,month,day,hour,min,rsec)
				  sec=int(rsec)
c				  write(time,'(i4,5i2.2)') year,month,day,hour,min,sec
				  write(time,'(i4,4i2.2,f6.3)') year,month,day,hour,min,rsec
				  write(write1,*) time,' ndet = ',ndet
				  write(write3,'(a,a,a,f6.1,a)')
     &             'wavetool ',base_out, ' -start '
     &             // time 
     &             // ' -duration ',extduration,' -wav_out_file SEISAN'
C WCC 3/2010: output picked stations and times
				  ncorrdet=ncorrdet+1
				  write (write4,'(a,i3,a)')
     &  '------------ ', ndet, ' stations ------------'
				  do k=i,j-1
				  	  if (flag(k).eq.1) then
					      hmsec=msec(k)
					      call sectim(hmsec,year,
     &                        doy,month,day,hour,min,rsec)
					      sec=int(rsec)
					      write (write4,'(a,1x,'
     &                      // 'i4,''/'',i2.2,''/'',i2.2,'' '','
     &                      // 'i2.2,'':'',i2.2,'':'',i2.2,'
     &                      // '1x,f6.1,1x,f6.1)') 
     &                           stat(k),year,month,day,hour,min,sec,
     &                           coda(k),ratio(k)
                       endif
				  enddo
C end WCC 3/2010
			  endif
          endif
      enddo
      
C		WCC added 3/2010
	  write(*,'(i8,a,f4.1,a,i3,a)') ncorrdet,
     &        ' events detected within ',
     &        netwindow, 's on ' , 
     &        netmindet , '+ stations'
      return
      end

C HEAPSORT ALGORITHM, FROM NUMERICAL RECIPES
C	Sorts an array ra(1:n) into ascending numerical order using the Heapsort algorithm. n is 
C	input; ra is replaced on output by its sorted rearrangement. 
C   MODIFIED FOR CONDET: ADD accompanying ARRAYS cas cac raa and rab,
C   and make ra double precision
	SUBROUTINE xhpsort(n,ra,cas,cac,raa,rab) 
	INTEGER n 
	double precision ra(n)
	real raa(n),rab(n)
	character*5 cas(n)
	character*4 cac(n)
	INTEGER i,ir,j,l 
	double precision rra
	real rraa,rrab
	character*5 ccas
	character*4 ccac
	integer iia
	
	if (n.lt.2) return 
C	The index l will be decremented from its initial value down to 1 during the
C   ÒhiringÓ (heap creation) phase. Once it reaches 1, the index ir will be
C   decremented from its initial value down to 1 during the
C   Òretirement-and-promotionÓ (heap selection) phase. 
	l=n/2+1 
	ir=n 
10 	continue 
	if (l.gt.1) then
		l=l-1 
		rra=ra(l) 
		ccas=cas(l)
		ccac=cac(l)
		rraa=raa(l)
		rrab=rab(l)
	else
		rra=ra(ir) 
		ccas=cas(ir)
		ccac=cac(ir)
		rraa=raa(ir)
		rrab=rab(ir)
		ra(ir)=ra(1)  
		cas(ir)=cas(1)
		cac(ir)=cac(1)
		raa(ir)=raa(1)
		rab(ir)=rab(1)
		ir=ir-1 
		if (ir.eq.1) then  
			ra(1)=rra 
			cas(1)=ccas
			cac(1)=ccac
			raa(1)=rraa
			rab(1)=rrab
			return 
		endif 
	endif 
	i=l
	j=l+l 
20	if(j.le.ir)then
		if (j.lt.ir) then 
		if (ra(j).lt.ra(j+1)) j=j+1 
		endif 
		if(rra.lt.ra(j)) then 
			ra(i)=ra(j) 
			cas(i)=cas(j)
			cac(i)=cac(j)
			raa(i)=raa(j)
			rab(i)=rab(j)
			i=j 
			j=j+j 
		else
			j=ir+1 
		endif 
		goto 20 
	endif 
	ra(i)=rra
	cas(i)=ccas
	cac(i)=ccac
	raa(i)=rraa
	rab(i)=rrab
	goto 10 
	END 

      subroutine netdet_old(read1,write1,write3,pre_event,extduration,
     &     netwindow,netmindet,base_out)

      implicit none
      character*5 base_out
      integer read1,write1,write3
      real pre_event,extduration
      real netwindow
      integer netmindet
      integer i,j,n,ind
      integer max
      parameter (max=1000000)
      double precision msec(max),hmsec
      logical used(max)
      character*5 stat(max),hstat
      integer year,month,day,hour,min,sec,doy
      real rsec
      integer ndet
      character*14 time

c
c read data    
c
      n=1
10    continue
      read(read1,'(a5,1x,4x,1x,i4,5i2)',end=50) 
     &   stat(n),year,month,day,hour,min,sec
      rsec=float(sec)
      call timsec(year,month,day,hour,min,rsec,msec(n))
      used(n)=.false.
      n=n+1
      if (n.ge.max) then
        write(*,*) ' dimensions exceeded, stop '
        stop
      endif
      goto 10

50    continue
      n=n-1
      write(*,*) ' number of detections: ',n
c
c sort data
c
      do i=1,n-1
        ind=i
        do j=i+1,n
          if (msec(j).lt.msec(ind)) then
            ind=j
          endif
        enddo
c
c swap
c
        if (ind.gt.0) then
          hstat=stat(i)
          hmsec=msec(i)
          stat(i)=stat(ind)
          msec(i)=msec(ind)
          stat(ind)=hstat
          msec(ind)=hmsec
        endif
      enddo

c
c now check if enough different stations in time window
c
      ndet=0
      do i=1,n
c       write(25,*) ' i ',i,netwindow
        j=i+1
        ndet=0
        if (.not.used(i)) then
         ndet=1
         if (i.lt.n) then
          do while(msec(j).le.msec(i)+netwindow.and.
     &             j.le.n) ! lot 02/10/2008
c         write(25,*) j,stat(j),msec(j),i,stat(i),msec(i)+netwindow,
c    &     msec(i)+netwindow-msec(j),ndet
            if (stat(j).ne.stat(i).and..not.used(j)) then
              ndet=ndet+1
c         write(25,*) ' j ',j,msec(j),msec(i)+netwindow,
c     &     msec(i)+netwindow-msec(j),ndet
              used(j)=.true.
              used(i)=.true.
            endif
            j=j+1
          enddo
         endif
        endif
c        write(*,*) ' ndet ',ndet
        if (ndet.ge.netmindet) then
          hmsec=msec(i)-pre_event
          call sectim(hmsec,year,
     &           doy,month,day,hour,min,rsec)
          sec=int(rsec)
          write(time,'(i4,5i2.2)') year,month,day,hour,min,sec
          write(*,*) time,' ndet = ',ndet
c          write(write3,*) 'wavetool -cwav -start '
                    write(write3,'(a,a,a,f6.1,a)')
c     &                'wavetool -cwav -start '
     &                'wavetool ',base_out,' -start '
     &                // time 
     &     // ' -duration ',extduration,' -wav_out_file SEISAN'
        endif
      enddo
      
      return
      end

