program spectrum
c********************************************************
c  Version 5-16-2010  M. J. Hinich
c  Periodogram, spectrum, & auto coherencies
c  File fname.out is always opened where 'fname' is the name of the data
c  This file contains the parameters used and the statistics of the data,
c  the residuals of an AR fit if used and the modulation signal.
c
c  If segment length=0 then the whole data file is one sgement & then
c   File fname.graph - a data structure file with the spectrum, signal coherence
c   spectrum, its probability value and the periodogram of the mean frame
c
c   File fname-signal.out - coherent part of the mean frame
c   File fname-modulation.out - modulation signal of the periodicity
c
c  If segment length>1 then
c   File fname.spectrogram - a data structure file with the spectra of each segment
c
c   File fname.cohprogram - a data structure file with the signal coherence spectra
c   and  probability values of each segment
c=================================================
c  Parameters in control file Runspectrum.cnl. Use spaces between entries!
c
c  A delimiter symbol * must be placed on each line.
c  Place comments on the lines of runspectrum.cnl file after the symbol *.
c  Runspect.cnl allows runs a number of control files with different file
c  names containing different data files and parameters.
c  The parameters for each run are echoed in filename.out where filename
c  is the name of the data file used.
c
c  The data file contains nc time series in column form where the data
c  format is set in the second line of the data file. If nc > 1 then each
c  column is read in as a separate run and the output files are numbered
c  from 1 to nc.
c
c===================================================================

c  Lines for the control file Runspectrum.cnl
c  The parameters & commands are found using key words whose first letter
c  may be either upper or lower case. The rest of the key word must be in
c  lower case. The symbol ! is a wildcard for the value to be set.
c  Examples: !nrun below is the wildcard for the number of runs of the program
c  using different cnl files called by runspectrum.cnl. The integer to the right
c  of the = sign in runspectrum.cnl is read by the program and sets the number
c  of runs.
c  In any *.cnl called by runspectrum.cnl, !r is the wildcard for the resolution
c  bandwidth. The real number to the right of !r= is the resolution bandwidth
c  that is used.
c
c  ALWAYS use 0.* rather than .* since the program searches for the number
c  to the left of the period in the parameter.
c
c  A delimiter symbol * must be placed on each line.
c  Use the keywords followed by an = sign and the setting. Some settings are
c  character strings while others are real numbers. A character string is
c  simply a sequence of characters that usually English words.
c
c  Comments on the lines of the cnl file are placed to the right of the
c  symbol *. These are not read by the program.
c
c===================================================================
c
c  Number of runs=!nrun   (!nrun - integer wildcard)
c
c  There has to be !nrun cnl files in runspectrum.cnl
c
c============================
c
c  Coherence threshold=!rth
c
c  !rth - threshold for computing the coherent part of the mean frame
c  The kth FFT bin is zeroed out if kth signal coherence < th.
c  The coherent part of the mean frame is in the file *.signal where * is the
c  data file name read from the control file read below.
c
c============================
c
c  These 2 entries are followed by the line:
c
c  # */ End of Runspect.cnl parameter reads
c
c  The program searches for the delimiter symbol # to end the reading of these
c  control parameters.
c
c  Below the delimiter line type the list of !nrun lines with the path & names
c  of the control files that pass parameters to the program.
c
c  For example if !nrun=2 & the control files are in the directory cnl:
c  \cnl\spectrum-fil1.cnl             * cnl file
c  \cnl\spectrum-fil2.cnl             * cnl file
c
c============================
c
c  The loop on runs k = 1,...,!nrun begins for !nrun control files.
c  The data file names for each control file & run parameters are now read.
c
c  The parameter lines are followed by the line
c  # * End of run parameter reads from the control file
c
c  The program searches for the delimiter symbol # to end the reading of the
c  parameters for each control file set in Runspect.cnl
c
c=========================
c
c  File name=data file name with path
c
c============================
c
c  Sampling unit=!unit
c
c  kHz or Hz to interpret the sampling rate as a multiple of kHz or Hz
c  You can also use the lower case letters k or h.
c
c  If khz is used the spectrum output unit is kHz, or if hz is used the
c  spectrum output unit is Hz   OR
c
c  a 10 character unit in lower case to interpret the sampling rate
c  as a multipled of this unit. The output will be in periods using that
c  unit as for example  unit='sec'.
c
c=======================================================================
c
c  Resolution bandwidth=!rb  ( !r - real wildcard)
c
c  This bandwidth is in Hz if the sampling unit is in frequency. The frame
c  length lb used to compute the spectrum is then inversely proportional to r.
c  If the sampling unit is NOT frequency then lb=rb.
c
c=========================
c  NOTE: sr, the 2nd number on line 2 of the data file is the sampling
c   rate in kHz or Hz, OR sr is a sampling interval with the unit determined by
c   the unit determined by the characters placed after unit=
c=======================================================================
c
c  Lower frequency=!rl  upper frequency=!ru
c
c  The unit is kHz or Hz depending on the frequency unit used
c
c  Longest period=!rl  shortest period=!ru if a unit is used
c
c  (!rl,!ru) is the bandwidth of the computed spectrum
c
c  If !rl=0 then the longest period is the resolution bandwidth period
c
c  If !ru=0 then the normalized upper frequencey is 0.5  or the
c  normalized shortest period is lb/2
c
c=========================
c
c  File read start=!ns  file read end=!ne
c
c  !ns - index of datum first read, !ne - index of datum last read
c  If !ne = 0, all the data is read.
c
c=========================
c
c  If you want to prewhiten with an AR(np) fit type the character string on
c  a line in the cnl file:
c   prewhiten with AR(!np) fit
c
c  !np - starting number of parameters of a recursive least squares AR
c  fit if np > 0. If np = 0 then this step is bypassed.
c
c  No AR fit is made if this character string is absent in the lines above
c  delimiter line that starts with the symbol # or if AR(0) is set on the line.
c
c  The residuals of the AR(np) are white if np is sufficiently large.
c  Using an AR(np) fit to generate white residuals is a linear operation
c  on the raw data.
c
c  The AR model is estimated in the subroutine far that outputs the
c  po < np used as the end the iterations.
c
c============================
c
c  If you want to detrend the data by fitting a curvilinear trend using the
c  first four Legendre polynomials type the following word on a line in the
c  cnl file:
c   detrend
c
c=========================
c
c  Prob threshold=!rpv
c
c  !rpv - probability level threshold for the t statistics for the AR
c  coefficients of each AR fit in the iterations or each OLS fit in the
c  Legendre polynomial detrending fit.

c  If a coefficient has a t value whose absolute value has a probability
c  Pr(abs(t)) < 1 - rpv, it is deleted in the next OLS fit.
c  The covariance matrix is adjusted for the subset of lagged values used.
c
c  The AR model is estimated in the subroutine far which outputs the
c  po < np used as the end of at most three iterations.
c
c  If rpv = 1 then the full AR(np) model is estimated.
c
c============================
c
c  Spectrogram segment length=!lw
c
c  The sampled signal is divided into overlapping segments of length !lw and
c  the power spectrum & signal coherence spectrum is computed for each segment.
c  The stacked log spectra are written to a file fname-spectrogram.out where
c  fname is the name of the data file and the stacked signal coherencies
c  probabilities are written to a file fname-cohprobgram.out.
c  The frequency or period is written on the last line of the data structure.
c
c  The same segment length is used for all the runs.
c
c  NOTE: lw must be < or = N/3 where N is the number of data rows in the data
c  file.
c
c  If length=0 then these spectrograms are not computed and the whole sample
c  is used once in computing the statistics. The modulation signals for each
c  column of data is written to a file fname-modulation.out.
c
c============================
c
c  Overlap=!ov
c
c  Each segmenet is overlapped by ov data points. If !ov=0 the segments are
c  not overlapped. If the segment length=0 then ov=0.
c
c============================
c
c  % taper=!rpt
c
c  !rpt = % taper of frame for sidelobe reduction (0. < rpt < 25.)
c  If !rpt = 0 or this line is omitted then the frames are NOT tapered.
c
c============================
c
c   Spectral smoothing bandwidth=!smb
c
c    The frame averaged spectrum is smoothed by a truncated cosine taper
c    whose bandwidth is !smb in units of Hz/milsec.
c    If the character string 'smoothing bandwidth' is not typed on a line
c    in the control file OR if smb = 0. then the spectrum is NOT smoothed.
c
c    smb MUST be > resolution bandwidth rb. If smb < rb, then smb = 0.
c
c=========================
c
c  Frame padding=!npad
c
c   !npad - number of extended frames padded by zeros to allow fine tuning
c   of the resolution bandwidth. For example if npad =1 & the frame length
c   is 100 then the frame is padded by 100 zeros to make a frame of length
c   of 200. The units are adjusted so that the output has correct units.
c   If !npd=0 or this linw ia omitted then no padding is done.
c
c=========================
c
c  Data on File 2
c  Data file format for the first two lines (records)
c  Line 1: A character string with the data file name and code (a80)
c  Line 2: rows=!nr columns=!nc  sampling rate (or interval)=!sr
c   format=(format)  unit=!
c   !nr - observations, !nc - data columns sr - sampling rate in kHz or Hz or
c   multiples of the unit set in the cnl file.
c   Data format is a character*50 string including ( )
c   unit=kHz or khz if that is the sr unit
c   unit=Hz or hz if that is the sr unit  OR unit=unit in time for sr
c
c   Example string: (4(f7.2,2x),2x,a7)
c   The time mark must be less than 26 characters long.
c   It is the last column of the data file.
c
c   Example string: (a25,1x,5(f7.2,1x))
c   The time mark must be less than 26 characters long.
c   It is the first column of the data file.
c
c===========================================
      integer run,nr,lw,ov
      real gs,rt
      character*700 name
      character*80 buf,cnl
      character*50 charead
      character*20 par
      character*10 su
      character*2 delim
      logical exs
c**********************
      open(109,file='Spectrum-error.out')
c  Inquire if file exists
      inquire(file='Runspectrum.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/'' Runspect.cnl file does not exist'')')
       write(109,'(/'' Runspectrum.cnl file does not exist'')')
       stop
      endif
c  Open Runspect.cnl
      open(1,file='Runspectrum.cnl',err=1,iostat=io1,status='old')
c  Start time
      call clock(gs,0,3)
c-------
      ip=0
      kt=0
      ia=0
      dowhile(kt<=17 .and. ia==0)
       kt=kt+1
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)',end=2,err=3,iostat=io1) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in Runspect.cn
     *l'
        write(109,*)'Place comment marker * after each line in Runspect.cn
     *l'
        write(6,'(a70)')buf
        write(109,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  No. of runs
      par='uns'
      delim(1:2)='= '
      run=numread(par,name,delim,ia,ier,io)
c  Coherence threshold
      par='ce threshold'
      th=rnumread(par,name,delim,ia,ier,io)
c  Trap th
      if(ier>0) then
       write(6,*) 'Coherence threshold ',th,' is not in the cnl file'
       write(109,*) 'Coherence threshold ',th,' is not in the cnl file'
      endif
      if(th<0 .or. th>=1.) then
       write(6,*)'Coherence threshold ',th,' is not in unit interval'
       write(109,*)'Coherence threshold ',th,' is not in unit interval'
       stop
      endif
c  Multiple runs
      do nr=1,run
c  Read control file
       read(1,'(a80)',end=3,err=4,iostat=io2)cnl
       ig=index(cnl,'*')
       cnl=cnl(:ig-1)
       call top(cnl,th,gs,run,nr)
      enddo
c  End of loop on runs
      stop
  1   write(6,*)'**** Error on opening Runspect.cnl ',io1
      stop
  2   write(6,*)'End of file on runspectrum.cnl reads'
      write(109,*)'End of file on runspectrum.cnl reads'
  3   write(6,*)'Place the symbol # after the 4 top reads ',io1
      write(109,*)'Place the symbol # after the 4 top reads ',io1
      stop
  4   write(6,*)'**** Error on reading name of cnl file ',nr,io2
      write(6,*)cnl
      write(109,*)'**** Error on reading name of cnl file ',nr,io2
      write(109,*)cnl
c-------
      stop
      end
c
c===========================================
c
      subroutine top(cnl,th,gs,run,nr)
c********************************************************
c  Runs each control file in the master Runspect.cnl
c========================================================
      parameter(mt=5,ndy=16,ny=2010)
      integer run,p,lw,ov,seg,mt,ndy,ny
      real gs,rt
      character*700 name
      character*120 fn
      character*80 buf,cnl
      character*50 charead,fmt,fname
      character*20 par
      character*10 su
      character*2 delim
      character dtnd,mk
      logical exs
      intent(in) cnl,th,gs,nr,run
c**********************
      inquire(file=cnl,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'cnl file does not exist ',cnl
       write(109,'(/)')
       write(109,*)'cnl file does not exist ',cnl
       return
      endif
c  Open control file
      open(11,file=cnl,err=1,iostat=io1,status='old')
c--------
      ip=0
      kt=0
      ia=0
      dowhile(kt<=20.and.ia==0)
       kt=kt+1
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(11,'(a80)',end=6)buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line in the cnl fil
     *e'
        write(6,'(a70)')buf
        write(109,*)'Place comment marker * after each line in the cnl fil
     *e'
        write(109,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
c  End of reads
      enddo
      delim(1:2)='='
c  Data file name from control file
      par='ile name'
      fn=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,*) 'Error in read of file name in cnl file'
       write(109,*)'Error in read of file name in cnl file'
       stop
      endif
      inm=index(fn,'.')
      if(inm==0) then
       write(109,*)'Data file ',fn,' has no extension'
       stop
      endif
c  Inquire if file exists
      inquire(file=fn,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'Input file ',fn(:inm+5),' does not exist'
       write(109,'(/)')
       write(109,*)'Input file ',fn(:inm+5),' does not exist'
       stop
      endif
c-------
c  Open data file
      open(2,file=fn,err=3,iostat=io3,status='old')
c  Parse for up to 5 path switches
      ib=index(fn,'\')
      if(ib>0) then
       ic=index(fn(ib:),'\')
       ib=ib+ic
       if(ic>0) then
        ic=index(fn(ib:),'\')
        ib=ib+ic
        if(ic>0) then
         ic=index(fn(ib:),'\')
         ib=ib+ic
         if(ic>0) then
          ic=index(fn(ib:),'\')
          ib=ib+ic
         endif
        endif
       endif
        else
         ib=1
      endif
      ia=index(fn,'.')
      fname=fn(ib:ia)
      nfn=ia-ib+1
      fname(nfn:nfn+10)='-Spect.out'
c  Open spectrum output file
      open(3,file=fname(:nfn+10),err=2,iostat=io2)
c  Name of program
      write(3,'(11x,''Signal Coherence and Spectrum'')')
      write(3,'(12x,''Melvin Hinich''/)')
c  Sampling unit
      par='unit'
      buf=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,'(/'' Error in read of sampling unit '',a15)')buf
       write(109,'(/'' Error in read of sampling unit '',a15)')buf
       stop
      endif
      su=buf(:15)
      ia=index(buf,'kHz')
      ib=index(buf,'khz')
      ic=index(buf,'Hz')
      ij=index(buf,'hz')
      if(ia>0.or.ib>0) then
c  kHz
       iunit=1
      endif
      if((ic>0.or.ij>0).and.(ia==0.and.ib==0)) then
c  Hz
       iunit=2
      endif
      if(ia==0.and.ib==0.and.ic==0.and.ij==0) then
c  Unit set in cnl file
       iunit=0
      endif
c-------
c  Read data file header
      call readhead(buf,nfl,nc,sr,fmt,iunit,2,3)
      write(3,'(1x,a79)')buf
c-------
      if(iunit>0) then
c  Lower frequency
       par='wer freq'
       fl=rnumread(par,name,delim,ia,ier,6)
       if(ier==1) then
        write(3,'(/'' Type lower freq rather than longest period in sett
     *ing the band in %.cnl'')')
        write(6,'(/'' Type lower freq rather than longest period in sett
     *ing the band in %.cnl'')')
        write(109,'(/'' Type lower freq rather than longest period in sett
     *ing the band in %.cnl'')')
        stop
       endif
c  Upper frequency
       par='per freq'
       fu=rnumread(par,name,delim,ia,ier,6)
       if(ier==1) then
        write(3,'(/'' Type upper freq rather than shortest period in set
     *ting the band in %.cnl'')')
        write(6,'(/'' Type upper freq rather than shortest period in set
     *ting the band in %.cnl'')')
        write(109,'(/'' Type upper freq rather than shortest period in set
     *ting the band in %.cnl'')')
        stop
       endif
        elseif(iunit==0) then
c  Longest period
         par='gest period'
         fl=rnumread(par,name,delim,ia,ier,6)
         if(ier==1) then
          write(3,'(/'' Type longest period rather than lower freq in se
     *tting the band in the cnl file'')')
          write(6,'(/'' Type longest period rather than lower freq in se
     *tting the band in the cnl file'')')
          write(109,'(/'' Type longest period rather than lower freq in se
     *tting the band in the cnl file'')')
          stop
         endif
c  Shortest period
         par='test period'
         fu=rnumread(par,name,delim,ia,ier,6)
         if(ier==1) then
          write(3,'(/'' Type shortest period rather than upper freq in s
     *etting the band in %.cnl'')')
          write(6,'(/'' Type shortest period rather than upper freq in s
     *etting the band in %.cnl'')')
          write(109,'(/'' Type shortest period rather than upper freq in s
     *etting the band in %.cnl'')')
          stop
         endif
      endif
c-------
c  Resolution bandwidth
      par='ion bandwidth'
      rb=rnumread(par,name,delim,ia,ier,6)
      if(ier==4) then
       write(3,*)'Place a period in the resolution bandwidth value'
       write(6,*)'Place a period in the resolution bandwidth value'
       write(109,*)'Place a period in the resolution bandwidth value'
       ia=index(name,'idth')
       write(3,*)name(ia+5:ia+15)
       write(6,*)name(ia+5:ia+15)
       write(109,*)name(ia+5:ia+15)
       stop
      endif
      if(rb<=0.0) then
       write(3,*)'Resolution bandwidth < 0'
       write(6,*)'Resolution bandwidth < 0'
       write(109,*)'Resolution bandwidth < 0'
       stop
      endif
c-------
      dtnd='n'
c  Curvilinear detrending
      ia=index(name,'trend')
      if(ia>0) then
       dtnd='y'
      endif
c-------
      p=0
      ia=index(name,'rewhite')
c  Prewhiten by an AR(p) fit
      if(ia>0) then
       ib=index(name(ia:ia+20),'AR')
       ic=index(name(ia:ia+20),'ar')
       if(ib>0) then
        par='R'
        delim(1:2)='()'
        p=numread(par,name(ib:),delim,ia,ier,6)
         elseif(ic>0) then
          par='r'
          delim(1:2)='()'
          p=numread(par,name(ic:),delim,ia,ier,6)
       endif
      endif
      delim(1:2)='= '
      if(p>0.or.dtnd=='y') then
c  Probability threshold for AR t values
       par='ty threshold'
       pv=rnumread(par,name,delim,ia,ier,6)
       if(ier>0) then
        write(6,*)'Error in read of prob threshold in cnl file'
        write(109,*)'Error in read of prob threshold in cnl file'
        stop
       endif
c  Trap pv
       if(pv<=0.) then
        write(3,*)'Prob level pv ',pv,' < 0'
        write(3,*)'Set pv = 0 if you do not want to prewhiten'
        write(6,*)'Prob level pv ',pv,' < 0'
        write(6,*)'Set pv = 0 if you do not want to prewhiten'
        write(109,*)'Prob level pv ',pv,' < 0'
        write(109,*)'Set pv = 0 if you do not want to prewhiten'
        stop
       endif
       if(pv>1) then
        write(3,*)'Prob level pv ',pv,' > 1'
        write(6,*)'Prob level pv ',pv,' > 1'
        write(109,*)'Prob level pv ',pv,' > 1'
        stop
       endif
      endif
c-------
c  % taper
      pt=0.0
      par='taper'
      pt=rnumread(par,name,delim,ia,ier,6)
      if(ier>0) then
       pt=0.0
      endif
c-------
c  Spectrum smoothing bandwidth
      smb=0.0
      par='hing bandwidth'
      smb=rnumread(par,name,delim,ia,ier,6)
      if(ier==4) then
       write(3,*)'Place a period in the spectrum smoothing bandwidth val
     *ue'
       write(6,*)'Place a period in the spectrum smoothing bandwidth val
     *ue'
       write(109,*)'Place a period in the spectrum smoothing bandwidth val
     *ue'
       ia=index(name,'idth')
       write(3,*)name(ia+5:ia+15)
       write(6,*)name(ia+5:ia+15)
       stop
        elseif(ier==0) then
         par='hing Bandwidth'
         smb=rnumread(par,name,delim,ia,ier,6)
         if(ier==4) then
          write(3,*)'Place a period in the spectrum smoothing bandwidth
     * value'
          write(6,*)'Place a period in the spectrum smoothing bandwidth
     * value'
          write(109,*)'Place a period in the spectrum smoothing bandwidth
     * value'
          ia=index(name,'idth')
          write(3,*)name(ia+5:ia+15)
          write(6,*)name(ia+5:ia+15)
          write(109,*)name(ia+5:ia+15)
          stop
         endif
      endif
      if(smb>rb.and.rb>0) then
c  No. of bins to smooth
       nw=nint(smb/rb)
       nw=2*(nw/2)+1
        else
         nw=0
      endif
c  Spectrogram segment length
      lw=0
      par='ength'
      lw=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,*)'Error in read of segment length in cnl file'
       write(109,*)'Error in read of segment length in cnl file'
       stop
      endif
c  Segment overlap
      ov=0
      par='rlap'
      ov=numread(par,name,delim,ia,ier,io)
      if(lw==0) then
c  No overlap
       ov=0
      endif
c  Interpolation padding
      npd=1
      par='adding'
      npd=numread(par,name,delim,ia,ier,6)
      if(ier==1) then
       npd=1
      endif
c-------
c  Trap npd
      if(npd<0) then
       write(3,*)'No. of pads ',npd,' < 0'
       write(6,*)'No. of pads ',npd,' < 0'
       write(109,*)'No. of pads ',npd,' < 0'
       stop
      endif
      if(npd==0) then
       npd=1
      endif
c-------
      if(iunit==0) then
c  Trap short period<sr
       if(fu<2.0*sr.and.fu>0.0) then
        write(6,'(/2x,''Shortest period '',g12.6,'' < 2*sampling interva
     *l '',g12.6/)')fu,2.0*sr
        write(109,'(/2x,''Shortest period '',g12.6,'' < 2*sampling inter
     *val '',g12.6/)')fu,2.0*sr
        stop
       endif
c  Trap fl=0 for low band period
       if(fl<=0.0) then
        fl=1.e6
       endif
c  Invert period to frequency
       fl=1.0/fl
       if(fu>0.0) then
        fu=1.0/fu
         else
          fu=0.5
       endif
      endif
      if(iunit>0) then
       if(fu>sr/2.0.and.fu>0.0) then
        write(6,'(/2x,''Upper frequency '',f12.3,'' > sampling rate/2 ''
     *f12.3/)')fu,sr/2.0
        write(109,'(/2x,''Upper frequency '',f12.3,'' > sampling rate/2
     * ''f12.3/)')fu,sr/2.0
        stop
       endif
       if(fu==0) then
        fu=0.5*sr
       endif
      endif
      if(fu>0.0.and.fl>fu) then
       write(3,*)'Upper band in %.cnl ',fu,' < lower band ',fl
       write(6,*)'Upper band in %.cnl ',fu,' < lower band ',fl
       write(109,*)'Upper bw in %.cnl ',fu,' < lower band ',fl
       if(iunit==1) then
        write(3,*)'Check to see the (fl,fu) units are kHz'
        write(6,*)'Check to see the (fl,fu) units are kHz'
        write(109,*)'Check to see the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(3,*)'Check to see the (fl,fu) units are Hz'
          write(6,*)'Check to see the (fl,fu) units are Hz'
          write(109,*)'Check to see the (fl,fu) units are Hz'
           else
            write(3,*)'Check to see the (fl,fu) units are ',su
            write(6,*)'Check to see the (fl,fu) units are ',su
            write(109,*)'Check to see the (fl,fu) units are ',su
       endif
       stop
      endif
c  Trap taper %
      if(pt<0.or.pt>25.) then
       write(3,*)'Taper % not in range (0,25) ',pt
       write(6,*)'Taper % not in range (0,25) ',pt
       write(109,*)'Taper % not in range (0,25) ',pt
       stop
      endif
c-------
c  File read start
      is=1
      par='d start'
      is=numread(par,name,delim,ia,ier,6)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of start datum'
       write(3,*)par
       write(6,*)'Error ',ier,' in read of start datum'
       write(6,*)par
       write(109,*)'Error ',ier,' in read of start datum'
       write(109,*)par
       stop
      endif
c  File read end
      ie=0
      par='d end'
      ie=numread(par,name,delim,ia,ier,6)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of end datum'
       write(3,*) par
       write(6,*)'Error ',ier,' in read of end datum'
       write(6,*)par
       write(109,*)'Error ',ier,' in read of end datum'
       write(109,*)par
       stop
      endif
      if(is<1) then
       write(3,*)'IS in %.cnl is < 1 ',is
       write(6,*)'IS in %.cnl is < 1 ',is
       write(109,*)'IS in %.cnl is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(3,*)'IE < IS in %.cnl ',is,ie
         write(6,*)'IE < IS in %.cnl ',is,ie
         write(109,*)'IE < IS in %.cnl ',is,ie
          stop
      endif
c------
c  Trap nc
      if(nc<1) then
       write(3,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(6,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(109,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       stop
      endif
c  Trap file length
      if(ie>nfl) then
       write(3,*)'IE ',ie, ' > N in data file',nfl
       write(6,*)'IE ',ie, ' > N in data file',nfl
       write(109,*)'IE ',ie, ' > N in data file',nfl
       stop
      endif
      if(ie==0) then
       ie=nfl
      endif
c  True sample size
      n=ie-is+1
c-------
c  Increase sampling rate if data is interpolated
      psr=dfloat(npd)*dble(sr)
c-------
c  Frame length lb
      if(iunit==1) then
       lb=nint(sngl(dble(psr)*1.d3/dble(rb)))
      endif
      if(iunit==2) then
       lb=nint(sngl(dble(psr)/dble(rb)))
       psr=psr*1.e-3
c  Convert frequencies from Hz to kHz
       fl=dble(fl)*1.d-3
       fu=dble(fu)*1.d-3
      endif
      if(iunit==0) then
c  Time unit & decrease sampling interval if data is interpolated
       psr=dble(sr)/dfloat(npd)
       lb=nint(sngl(dble(rb)/dble(psr)))
      endif
c  Trap lb
      if(lb<1) then
       write(3,*)'Frame length = ',lb,' < 1'
       write(6,*)'Frame length = ',lb,' < 1'
       write(109,*)'Frame length = ',lb,' < 1'
       if(iunit>0) then
        write(3,*)'Decrease resolution bandwidth = ',rb,' n = ',n
        write(6,*)'Decrease resolution bandwidth = ',rb,' n = ',n
        write(109,*)'Decrease resolution bandwidth = ',rb,' n = ',n
         else
          write(3,*)'Increase resolution period = ',rb,' n = ',n
          write(6,*)'Increase resolution period = ',rb,' n = ',n
          write(109,*)'Increase resolution period = ',rb,' n = ',n
       endif
       stop
      endif
c  Trap lw
      if(lw>0.and.lw>n/2) then
       write(3,'(/7x,''ERROR - segment length > N/2 = '',i6)')n/2
       write(3,'(/7x,''Segment length set = '',i6/)')lw
       write(6,'(/7x,''ERROR - segment length > N/2 = '',i6)')n/2
       write(6,'(/7x,''Segment length set = '',i6/)')lw
       write(109,'(/7x,''ERROR - segment length > N/2 = '',i6)')n/2
       write(109,'(/7x,''Segment length set = '',i6/)')lw
       stop
      endif
      if(lw>0.and.lw<lb) then
       write(3,'(/7x,''ERROR - segment length < frame = '',i6)')lb
       write(3,'(/7x,''Segment length set = '',i6/)')lw
       write(6,'(/7x,''ERROR - segment length < frame = '',i6)')lb
       write(6,'(/7x,''Segment length set = '',i6/)')lw
       write(109,'(/7x,''ERROR - segment length < frame = '',i6)')lb
       write(109,'(/7x,''Segment length set = '',i6/)')lw
       stop
      endif
c-------
c  Convert to normalized values (0,.5)
      if(iunit>0) then
       fl=fl/psr
       if(fu>0.) then
        fu=fu/psr
         else
          fu=0.5
       endif
        else
         fl=fl*psr
         fu=fu*psr
      endif
c-------
c  Trap top band
      if(fu>0.5) then
       write(3,*)'Normalized upper band = ',fu, ' > 0.5'
       write(6,*)'Normalized upper band = ',fu, ' > 0.5'
       write(109,*)'Normalized upper band = ',fu, ' > 0.5'
       if(iunit==1) then
        write(3,*)'Check to see if the (fl,fu) units are kHz'
        write(6,*)'Check to see if the (fl,fu) units are kHz'
        write(109,*)'Check to see if the (fl,fu) units are kHz'
         elseif(iunit==2) then
          write(3,*)'Check to see if the (fl,fu) units are Hz'
          write(6,*)'Check to see if the (fl,fu) units are Hz'
          write(109,*)'Check to see if the (fl,fu) units are Hz'
           else
            write(3,*)'Check to see if the (fl,fu) units are ',su
            write(6,*)'Check to see if the (fl,fu) units are ',su
            write(109,*)'Check to see if the (fl,fu) units are ',su
       endif
       write(3,*)' Sampling unit ',su,' Sampling rate/interval ',psr
       write(3,*)'iunit = ',iunit
       write(6,*)' Sampling unit ',su,' Sampling rate/interval ',psr
       write(6,*)'iunit = ',iunit
       write(109,*)' Sampling unit ',su,' Sampling rate/interval ',psr
       write(109,*)'iunit = ',iunit
       stop
      endif
c  Trap p
      if(p<0) then
       write(3,*)'AR parameter ',p,' < 0. p set = 0'
       write(6,*)'AR parameter ',p,' < 0. p set = 0'
       write(109,*)'AR parameter ',p,' < 0. p set = 0'
       p=0
        elseif(p>(n/4)) then
         write(3,*)'AR parameter ',p,' > n/4. p set = 0'
         write(6,*)'AR parameter ',p,' > n/4. p set = 0'
         write(109,*)'AR parameter ',p,' > n/4. p set = 0'
         p=0
      endif
c-------
c  Trap spectral smoothing width
      nsmooth=nint(0.75*float(lb))
      if(nw>nsmooth) then
       write(3,*)'Spectral smoothing bandwidth is too large'
       write(3,*)'No. of bins to be smoothed ',nw,' > (3/4)lb= ',nsmooth
       write(6,*)'Spectral smoothing bandwidth is too large'
       write(6,*)'No. of bins to be smoothed ',nw,' > (3/4)lb= ',nsmooth
       write(109,*)'Spectral smoothing bandwidth is too large'
       write(109,*)'No. of bins to be smoothed ',nw,' > (3/4)lb= ',nsmoo
     *th
       stop
      endif
c-------
      write(3,'(/2x,''First datum used '',i8,2x,''Last datum used'',i8,
     *2x,''Sample size = '',i8)')is,ie,n
c-------
c  Increase sample & segment length=nsg if npd>1
      nn=npd*n
      if(lw==0) then
       nsg=nn
        elseif(lw>0) then
         nsg=npd*lw
          else
           write(3,'(/4x,''ERROR - segment length is wrong''/)')
           write(6,'(/4x,''ERROR - segment length is wrong''/)')
           write(109,'(/4x,''ERROR - segment length is wrong''/)')
           stop
      endif
      np=nn/lb
c  Trap no. of frames
      if(np<2) then
       write(3,*)'No. frames ',np,' < 2',' lb ',lb,' pads ',npd
       write(6,*)'No. frames ',np,' < 2',' lb ',lb,' pads ',npd
       write(109,*)'No. frames ',np,' < 2',' lb ',lb,' pads ',npd
       if(iunit>0) then
        write(3,*)'Increase resolution bandwidth = ',rb,' n = ',nn
        write(6,*)'Increase resolution bandwidth = ',rb,' n = ',nn
        write(109,*)'Increase resolution bandwidth = ',rb,' n = ',nn
         else
          write(3,*)'Decrease resolution period = ',rb,' n = ',nn
          write(6,*)'Decrease resolution period = ',rb,' n = ',nn
          write(109,*)'Decrease resolution period = ',rb,' n = ',nn
       endif
       stop
      endif
c  Seg - number of segments with or without overlap
      if(ov>0) then
c  Segment moved ov at each iteration
       seg=(nn-nsg)/ov+1
       if(ov*(seg-1)+nsg+1>nn) then
        seg=seg-1
       endif
        elseif(ov==0) then
c  No overlap
         seg=nn/nsg
      endif
c  Trap np
      if(seg<1)then
       write(3,'(/7x,''ERROR - Segment length is too large''/)')
       write(3,'(7x,''Segment = '',i6,'' Sample = '',i7/)')nsg,nn
       write(6,'(/7x,''ERROR - Segment length is too large''/)')
       write(6,'(7x,''Segment = '',i6,'' Sample = '',i7/)')nsg,nn
       write(109,'(/7x,''ERROR - Segment length is too large''/)')
       write(109,'(7x,''Segment = '',i6,'' Sample = '',i7/)')nsg,nn
       return
      endif
      if(seg==1) then
c  Open fname.signal for coherent part of the mean frame
       fname(nfn:nfn+7)='.signal'
       open(101,file=fname(:nfn+7),err=4,iostat=io4)
c  Open fname-modulation.out for modulations signals of the data columns
       fname(nfn:nfn+15)='-modulation.out'
       open(102,file=fname(:nfn+15),err=5,iostat=io5)
      endif
c-------
      call array(n,nw,ie,is,lb,nsg,npd,fl,fu,rb,psr,pt,p,fmt,su,fn,pv,th
     *,nc,iunit,ov,seg,dtnd,3)
      if(nc>0) then
       write(3,'(7x,''Number of Data Columns = '',i2/)')nc
      endif
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      return
  1   write(6,*)'**** Error on opening cnl file no. ',nr,io1
      write(6,*) cnl
      write(109,*)'**** Error on opening cnl file no. ',nr,io1
      write(109,*)cnl
      stop
  2   write(6,*)'**** Error on opening output file ',io2
      write(109,*)'**** Error on opening output file ',io2
      stop
  3   write(6,*)'**** Error on opening data file ',io3
      write(6,*)fn
      write(109,*)'**** Error on opening data file ',io3
      write(109,*)fn
      stop
  4   write(6,*)'**** Error on opening coherent part .sig file ',io4
      write(109,*)'**** Error on opening coherent part .sig file ',io4
      stop
  5   write(6,*)'**** Error on opening modulations file ',io5
      write(109,*)'**** Error on opening modulations file ',io5
      stop
  6   write(6,'(7x,''ERROR''/)')
      write(6,*)'End of file in the read of the cnl parameters'
      write(6,*)'Place the symbol # at the end of the parameter list in
     * the cnl file'
      write(6,'(a80)')cnl
      write(109,*)'End of file in the read of the cnl parameters'
      write(109,*)'Place the symbol # at the end of the parameter list i
     *n the cnl file'
      write(109,'(a80)')cnl
      stop
      end
c
c===========================================
c
      subroutine array(n,nw,ie,is,lb,nsg,npd,fl,fu,rb,sr,pt,p,fmt,su,fn,
     *pv,th,nc,iunit,ov,seg,dtnd,io)
c********************************************************
c  nsg - segment lenth, seg - no. of segments, ov - overlap
c======================
      allocatable::x,gr,sp,mod,a,v,w,y,z,ao,ik,f,g,rt,ts
      real x(:),y(:),z(:),a(:),gr(:),sp(:),mod(:),w(:),fl,fu
      real*8 v(:)
      complex ao(:),rt(:)
      complex*16 f(:),g(:)
      integer ik(:),p,n,is,ie,lb,npd,lw,ov,seg
      character*120 fn
      character*80 buf
      character*50 fmt
      character*25 ts(:)
      character*10 su
      character dtnd,mk
      intent(in) n,nw,ie,is,lb,nsg,fl,fu,sr,npd,pt,fmt,su,p,pv,th,nc,iun
     *it,ov,rb,seg,dtnd,io
c**********************
c  Set il - bin for fu
      dlb=dfloat(lb)
      il=nint(sngl(dble(fl)*dlb))
      il=max(il,1)
      n2=lb/2
c  Set iu - upper band integer limit
      iu=nint(sngl(dble(fu)*dlb))
      iu=min(iu,n2)
      if(il>=iu) then
       write(3,'(2x,''Lower bin no. '',i4,2x,''Upper bin no. '',i6/)')il
     *,iu
       write(6,'(2x,''Lower bin no. '',i4,2x,''Upper bin no. '',i6/)')il
     *,iu
       write(109,'(2x,''Lower bin no. '',i4,2x,''Upper bin no. '',i6/)')
     *il,iu
       if(iunit>0) then
        write(3,'(2x,''Decrease resolution bandwidth '',f9.4)')rb
        write(3,'(/4x,''Frame length = '',i7)')lb
         else
          write(3,'(2x,''Increase resolution bandwidth '',f9.4)')rb
          write(3,'(/4x,''Frame length = '',i7)')lb
       endif
       if(iunit>0) then
        write(6,'(2x,''Decrease resolution bandwidth '',f9.4)')rb
        write(6,'(/4x,''Frame length = '',i7)')lb
         else
          write(6,'(2x,''Increase resolution bandwidth '',f9.4)')rb
          write(6,'(/4x,''Frame length = '',i7)')lb
       endif
       if(iunit>0) then
        write(109,'(2x,''Decrease resolution bandwidth '',f9.4)')rb
        write(109,'(/4x,''Frame length = '',i7)')lb
         else
          write(109,'(2x,''Increase resolution bandwidth '',f9.4)')rb
          write(109,'(/4x,''Frame length = '',i7)')lb
       endif
       stop
        elseif(il>=iu .and. iu==n2) then
         write(io,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
         write(6,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
        write(109,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
        stop
      endif
c  No. of frequency bins in band
      jf=iu-il+1
c-------
      nn=npd*n
      na=max(nn,nsg)
      n1=n*nc-1
      n2=na-1
      n3=nc*lb-1
      n4=lb-1
      n5=p-1
      n6=nsg-1
      n7=2*jf*seg-1
c-------
c  Allocate space
      allocate(x(0:n1),mod(0:n1),y(0:n2),z(0:n3),a(0:n4),sp(0:n4),w(0:n2
     *),ao(0:p),rt(0:p-1),ik(0:n5),f(0:n6),g(0:n4),gr(0:n7),v(0:4*nsg+14
     *),stat=ibad)
      if(ibad/=0) then
       write(3,*) 'Unable to allocate work space for arrays in array'
       stop
      endif
c-------
c  Parse format to see if there is an 'a' for time values
c  If so then set mk='m' to flag times
      mk='a'
      ia=index(fmt,'a')
      if(ia>0) then
       mk='m'
      endif
      if(mk=='m') then
       nmk=ie
       allocate(ts(0:nmk-1),stat=ibad)
       if(ibad/=0) then
        write(3,*) 'Unable to allocate work space for array ts'
        write(6,*) 'Unable to allocate work space for array ts'
        write(109,*) 'Unable to allocate work space for array ts'
        stop
       endif
        else
         nmk=1
         allocate(ts(0:0))
      endif
c  Write header for summary stats in file 3
      call wrh(su,sr,rb,n,ov,nsg,seg,npd,lb,pt,p,pv,th,nw,fl,fu,iunit,io
     *)
      call roc(x,sp,mod,a,w,y,z,ts,ao,rt,f,g,gr,v,ik,n,nw,lb,npd,fl,fu,i
     *s,ie,il,jf,pt,sr,su,nmk,p,mk,fmt,pv,rb,th,fn,na,nc,iunit,ov,nsg,se
     *g,dtnd,io)
c-------
      deallocate(x,sp,mod,a,w,y,z,ao,ik,f,g,gr,rt,ts,v)
      return
      end
c
c===========================================
c
      subroutine roc(x,sp,mod,a,w,y,z,ts,ao,rt,f,g,gr,v,ik,n,nw,lb,npd,f
     *l,fu,is,ie,il,jf,pt,sr,su,nmk,p,mk,fmt,pv,rb,th,fn,na,nc,iunit,ov,
     *nsg,seg,dtnd,io)
c********************************************************
c  Read data & loop on nc columns
c========================================================
      real x(n,nc),mod(n,nc),y(na),z(lb,nc),a(lb),gr(2*jf,seg),sp(lb),w(
     *na),fl,fu,spmax,spmin,spx,spn,amin,amax
      real*8 v(4*nsg+15),c(4)
      integer ik(p),in(4),nc,p,po,ov,seg,t,iunit,nsg,is,ie,jf,nw,lb,io
      complex ao(p+1),rt(p)
      complex*16 f(0:nsg-1),g(0:lb-1)
      character*120 fn
      character*50 fmt,fname
      character*40 id
      character*25 ts(nmk)
      character*10 su
      character*2 chno,co
      character dtnd,mk
      intent(in) n,nw,lb,fl,fu,is,ie,il,jf,npd,pt,sr,su,nmk,p,mk,fmt,pv,
     *rb,th,fn,na,nc,iunit,ov,nsg,seg,dtnd,io
c**********************
c  Read data
      call datc(x,ts,fmt,mk,nmk,is,ie,n,nc,2,io)
c  nsg - (padded) segment length
      if(ov==0) then
       novlap=nsg
        elseif(ov>0) then
c  Set pointer for overlap
         novlap=ov
          else
           write(6,*)'Overlap parameter ',ov,' is incorrect'
           write(3,*)'Overlap parameter ',ov,' is incorrect'
           stop
      endif
c-------
      ia=index(fn,'.')
c  Parse for up to 4 path switches
      ib=index(fn,'\')
      if(ib>0) then
       ic=index(fn(ib:),'\')
       ib=ib+ic
       if(ic>0) then
        ic=index(fn(ib:),'\')
        ib=ib+ic
        if(ic>0) then
         ic=index(fn(ib:),'\')
         ib=ib+ic
        endif
       endif
        else
         ib=1
      endif
      fname=fn(ib:ia)
      nfn=ia-ib+1
      n1=nc+6
      nc2=nc+nc+6
c  Max and min of spectrum for initial min & max
      amax=-1.0
      amin=1.e7
c  Loop on nc columns
      do kr=1,nc
c  Column kr name
       co=chno(kr)
       fname(nfn:nfn)='-'
       fname(nfn+1:nfn+2)=co
       rewind(2)
       read(2,'(a40)')id
       write(io,'(16x,12(''=''))')
       write(io,'(17x,''Column '',i3)')kr
       write(io,'(16x,12(''='')/)')
       if(seg==1) then
c  Open fname-co.graph
        close(n1+kr)
        fname(nfn+2:nfn+8)='.graph'
        open(n1+kr,file=fname(:nfn+8),err=1,iostat=io1)
        write(n1+kr,'(a40)')id
       endif
       if(seg>1) then
        close(n1+kr)
        fname(nfn+2:nfn+14)='.spectrogram'
        open(n1+kr,file=fname(:nfn+14),err=2,iostat=io2)
        write(n1+kr,'(a40)')id
        close(nc2+kr)
        fname(nfn+2:nfn+14)='.cohprobgram'
        open(nc2+kr,file=fname(:nfn+14),err=3,iostat=io3)
        write(nc2+kr,'(a40)')id
       endif
c-------
c  Shift kr column to y
       do t=1,n
        y(t)=x(t,kr)
       enddo
c  Compute statistics for input data
       call stat(n,y,am,sd,c3,c4,c6,smax,smin,io)
       write(io,'(60(''*''))')
       write(io,'(7x,''Descriptive Statistics of the Data = '',i8)')n
       call wrstat(io,am,sd,c3,c4,c6,smax,smin)
c  Normalize data
       do t=1,n
        y(t)=(y(t)-am)/sd
       enddo
c-------
       po=0
       if(dtnd=='y') then
c  Detrend
        call detrend(y,w,c,in,n,pv,po,io)
c  Residuals w =>y
        do t=1,n
         y(t)=w(t)
        enddo
       endif
c-------
       po=0
       if(p>0) then
c  Least squares AR(p) - residuals in array w, coefficients in v
        call far(y,v,a,w,ik,r2,sig,n,p,po,pv,ier,io)
        if(ier>0) then
         write(io,*)'Error in subroutine far ',ier
         write(6,*)'Error in subroutine far ',ier
         return
        endif
c  Compute statistics of the residuals
        call stat(n,w,am,sd,c3,c4,c6,smax,smin,io)
c-------
c  White test with exponent exp
        exp=0.33
        call port(w,n,exp,ct,io)
        if(po>0) then
         if(iunit>0) then
          srt=sr/float(npd)
           else
            srt=float(npd)*sr
         endif
         call outar(v,a,ao,rt,ik,r2,sig,su,srt,pv,po,iunit,io)
        endif
        write(io,'(60(''*''))')
        write(io,'(10x,''Descriptive Statistics of AR Residuals'')')
        call wrstat(io,am,sd,c3,c4,c6,smax,smin)
       write(io,'(/7x,''Whiteness test statistic p-value = '',f7.3/)')ct
c-------
c  Interpolate data for grid unit/npd if npd > 1 =>w
        if(npd>1) then
         call interp(w,y,f,v,lw,npd,io)
c  y =>w
         do t=1,npd*n
          w(t)=y(t)
         enddo
        endif
       endif
       if(p==0) then
        if(npd>1) then
         call interp(y,w,f,v,lw,npd,io)
          else
c  y =>w
           do t=1,n
            w(t)=y(t)
           enddo
        endif
       endif
c-------
c  Loop over segments
       do kp=1,seg
        spmax=0.0
        spmin=1.e7
c  Segment index
        km=(kp-1)*novlap
c  Shift kp segment to y
        do t=km+1,km+nsg
         y(t-km)=w(t)
        enddo
        call rn(y,gr,z,sp,mod,a,w,f,g,v,n,nw,lb,fl,fu,pt,sr,su,th,ty,kr,
     *na,nc,iunit,nsg,npd,seg,il,jf,kp,io)
        if(seg==0) then
         write(io,'(/7x,''p-value = '',f7.5)')1.0-ty
         call rmaxmin(sp,jf,spmax,spmin,amin,amax)
         write(io,'(/2x,''Maximum of Log Spectrum = '',f9.4,2x,''Minimum
     * = '',f9.4)')10.0*alog10(spmax),10.0*alog10(spmin)
        endif
        if(seg>1) then
         call rmaxmin(sp,jf,spx,spn,amin,amax)
         if(spx>spmax) then
          spmax=spx
         endif
         if(spn<spmin) then
          spmin=spn
         endif
c-------
c  Segment start time in minutes
         if(iunit==1) then
          start=float(km)/(60.0*sr)
           elseif(iunit==2) then
            start=float(km)/(60.0*sr*1.e3)
             elseif(iunit==0) then
              start=float(km)
         endif
        endif
        write(io,'(/'' Segment '',i5,2x,''Start time = '',f14.5,'' minut
     *es  p-value = '',f7.5)')kp,start,1.0-ty
       enddo
c  End of segments loop
       if(seg>1) then
c  Spectrogram & cohprobgram for column kr
        call wrtgrmc(gr,sr,su,lb,il,jf,seg,iunit,nc,kr,km,io)
        write(io,'(/2x,''Maximum of Log Spectrum = '',f9.4,2x,''Minimum
     *= '',f9.4)')10.0*alog10(spmax),10.0*alog10(spmin)
       endif
      enddo
c  End of column loop
      if(seg==1) then
c  Write coherent mean
       if(nc>1) then
        write(101,'(1x,''Coherent Part of the Mean Frame of the '',i2,'' C
     *olumns'')')nc
         else
          write(101,'(1x,''Coherent Part of the Mean Frame'')')
       endif
       do k=1,lb
        write(101,'(:,100(f17.7))')(z(k,kcm),kcm=1,nc)
       enddo
       call flush(101)
c  Output mod
       write(102,'('' Signal Modulations'')')
       nrws=(n/lb)*lb
       if(iunit==0) then
        write(102,'(''rows='',i7,'' columns='',i2,'' sampling interval=1
     *.0 format=(f14.5,1x) unit='',a10)')nrws,nc,su
       endif
       if(iunit==1) then
        write(102,'(''rows='',i7,'' columns='',i2,'' sampling rate='',f1
     *2.3,'' format=(f14.5,1x) unit=kHz'')')nrws,nc,sr
       endif
       if(iunit==2) then
        write(102,'(''rows='',i7,'' columns='',i2,'' sampling rate='',f1
     *2.3,'' format=(f14.5,1x) unit=Hz'')')nrws,nc,sr*1.e3
       endif
       do t=1,nrws
        write(102,'(:,10000(f14.5,1x))')(mod(t,ncm),ncm=1,nc)
       enddo
       call flush(102)
      endif
c-------
      write(io,'(60(''=''))')
      write(io,'(7x,''Header for Output in *.graph file''/)')
      if(iunit==1) then
       write(io,'(2x,''KHz, Log Spectrum (dB), Coherence, Chi-Square Pro
     *bability, Periodogram'')')
      endif
      if(iunit==2) then
       write(io,'(2x,''Hz, Log Spectrum (dB), Coherence, Chi-Square Prob
     *ability, Periodogram'')')
      endif
      if(iunit==0) then
       write(io,'(2x,''Period '',a10,'', Log Spectrum (dB), Coherence, C
     *hi-Square Probability, Periodogram'')')su
      endif
      write(io,'(/7x,'' Format: '',a50/)')fmt
c  Write large sample se
      bk=float(n/lb)
      cb=10.0/(alog(10.0)*sqrt(bk))
      write(io,'(7x,''Large Sample Standard Error = '',g9.4/)')cb
c-------
      call flush(io)
      return
  1   write(3,*)'**** Error opening spectrum.graph file ',io1
      stop
  2   write(6,*)'**** Error opening spectrogram file ',io2
      stop
      write(109,*)'**** Error opening spectrogram file ',io2
      stop
  3   write(6,*)'**** Error opening cohprobgram file ',io3
      stop
      write(109,*)'**** Error opening cohprobgram file ',io3
      stop
      end
c
c===========================================
c
      subroutine rn(y,gr,z,sp,mod,a,w,f,g,v,n,nw,lb,fl,fu,pt,sr,su,th,ty
     *,kr,na,nc,iunit,nsg,npd,seg,il,jf,kp,io)
c********************************************************
c  Compute stats =. graph for kr column, nsg - seg length, ty - stationarity cdf
c  Input: y - data array, su - unit, sr - sampling rate, lb - frame length
c  fl - lower freq, fu - upper freq, pt - % taper, th - coh thresh, kr - column
c  nw - no. of spectrum bins smoothed by convol, seg - segments, kp - seg index
c  gr(1,:) - sigcoh probs, gr(jf+1:) - dB spectrum for column kr
c  w - coherences, z - coherent part of mean kr frame, linear spectrum
c========================================================
      allocatable::b,r
      real y(na),gr(2*jf,seg),mod(n,nc),z(lb,nc),a(lb),sp(lb),w(na),b(:)
     *,r(:)
      real*8 v(4*nsg+15)
      integer seg
      complex*16 f(0:nsg-1),g(0:lb-1)
      character*10 su
      intent(in) y,n,nw,lb,fl,fu,sr,pt,su,th,kr,na,nc,iunit,nsg,npd,seg,
     *il,jf,kp,io
      intent(out) gr,mod,sp,ty,z
c**********************
      allocate(b(0:lb-1),r(0:n-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for b & r'
       write(6,*)'Unable to allocate work space for b & r'
       stop
      endif
c  a  - chi**2, b - |A|**2/lb, g - FFT mean frame, sp - spectrum, v - for lb
      call cohern(y,sp,mod,a,b,f,v,g,r,sg,fl,fu,nsg,lb,il,jf,nw,npd,pt,t
     *y,nc,kr,io)
      iu=il-1+jf
      call flush(io)
c  gr(1,:) - sigcoh probs, gr(jf+1:) - spectrum for column kr
c  w - coherences, z - coherent part of mean kr frame
      call grams(a,gr,sp,w,th,lb,il,jf,seg,nsg,kp,io)
      if(seg==1) then
c  z - coherent part of mean frame, sg - standard deviation of mean frame
       call cohsig(w,z,sg,f,g,v,lb,il,jf,sr,th,iunit,nc,seg,kr,io)
      endif
      call flush(io)
c  Output spectrum
      if(seg==1) then
       call wrtspec(gr,w,b,lb,il,jf,seg,su,sr,iunit,nc+6+kr,io)
      endif
c-------
      deallocate(b,r)
      return
      end
c
c===========================================
c
      subroutine cohern(x,sp,mod,a,b,f,v,g,r,sg,fl,fu,n,lb,il,jf,nw,npd,
     *pt,ty,nc,kr,io)
c********************************************************
c  a  - chi**2, b - |A|**2/lb, g - FFT of mean frame, sp - spectrum
c  v - initialized for lb, sg - standard deviation of mean frame
c  mod - modulation signals, kr - column, ty - stationarity test cdf
c========================================================
      real x(n),mod(n,nc),r(n),sp(lb),a(lb),b(lb),df,ty
      real*8 v(4*n+15)
      complex*16 f(0:n-1),g(0:lb-1)
      integer seg,t,n,lb,il,jf,nw,npd,kr,nc,n1,km1
      intent(in) x,n,lb,npd,fl,fu,il,jf,nw,pt,nc,kr
      intent(out) sp,mod,a,b,g,v,sg,ty
c**********************
      rl=float(lb)
      rp=float(npd)
      nrws=(n/lb)*lb
      scc=float(nrws)/(rl*rl)
c  a - squared cabs of FFT's of mean frame
      call amp(x,r,a,b,g,v,n,lb,il,jf,nt,pt,sg,sw,io)
c  Modulation r for column kr =>mod
      do t=1,n
       mod(t,kr)=r(t)
      enddo
c  Statistics of the modulation signal
       call stat(nrws,r,am,sd,c3,c4,c6,smax,smin,io)
       write(io,'(60(''*''))')
       write(io,'(7x,''Descriptive Statistics of the Modulation'')')
       call wrstat(io,am,sd,c3,c4,c6,smax,smin)
       write(io,'(60(''=''))')
c-------
c  Compute spectrum - b=w in spec
      call spec(r,b,sp,f,v,n,lb,npd,il,jf,nt,sw,io)
      chst=0.
      df=0.
c  Chi**2 loop
      do k=1,jf
c  Chi square
       ch=a(k)/sp(k)*rp
c  Add sinusoid's power back to spectrum
       sp(k)=sp(k)+a(k)/rl
c  |A|**2/lb =>b
       b(k)=a(k)/rl
c  Save chi**2 => a
       a(k)=ch
c  Accumulate chi squares
       chst=chst+ch
       df=df+1.0
      enddo
c  End of chi**2 looop
      chst=2.*scc*chst
      call cdchi(chst,2.0*df,ty,io)
c-------
c  Smooth spectrum if nw > 2
      if(nw>2) then
       ks=il-1
       iu=ks+jf
c  Shift spectrum values to sp(lb/2+1)
       do k=1,n2
        if(k<il) then
         sp(n2+k)=0.0
          elseif(k>=il .and. k<=iu) then
           sp(n2+k)=sp(k-ks)
            else
             sp(n2+k)=0.0
        endif
       enddo
       call convolv(sp(n2+1),ws,n2,nw,il,iu,ier,io)
       if(ier>0) then
        write(io,*) 'ier = ',ier,' from CONVOL'
       endif
c  Shift back smoothed spectrum into sp
       do k=1,jf
        sp(k)=sp(n2+ks+k)
       enddo
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine amp(x,r,a,b,g,v,n,lb,il,jf,nt,pt,sg,sw,io)
c********************************************************
c  Complex amplitudes of mean signal
c   il - bin no. of lower freq, jf - No. of frequencies used
c  Output: a - squared cabs of FFT's of mean frame
c   g - FFT of mean frame, v - initialized for lb, b - taper values
c   sg - standard deviation of mean frame, r - modulation signal
c========================================================
      real x(n),r(n),a(lb),b(lb)
      real*8 v(4*lb+15),z
      complex*16 g(0:lb-1)
      intent(in) x,n,lb,il,jf,pt
      intent(out) a,b,g,r,v,nt,sg,sw
c**********************
c  Average frame of x & subtract mean frame =>r
      call wave(x,r,a,n,lb,sg,io)
c  Initialize taper
      nt=0.
      sw=1.
      if(pt>0.) then
       call taper(b,lb,pt,sw,nt,io)
c  If no. of tapered values at ends < 1 set nt=0
       if(nt<1) then
        nt=0
        sw=1.0
       endif
      endif
c-------
c  Initialize v
      call cffti(lb,v)
c  Move a to complex g
      do j=1,lb
       g(j-1)=dcmplx(dble(a(j)))
      enddo
c-------
c  FFT of average frame
      call cfftf(lb,g,v)
c  Compute squared magnitude
      ks=il-1
      do k=1,jf
       kn=k+ks
c  Squared cabs(a(k))
       z=dreal(g(kn))*dreal(g(kn))+dimag(g(kn))*dimag(g(kn))
       a(k)=sngl(z)
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine wave(x,r,a,n,lb,sg,io)
c********************************************************
c  Subtracts mean frame a from x
c  Output: a - mean frame, r - signal minus mean frame, sg - standard dev.
c========================================================
      real x(n),r(n),a(lb)
      real*8 sm,sv
      intent(in) x,n,lb,io
      intent(out) a,r,sg
c*********************
      np=n/lb
c  Initialize signal array
      do k=1,lb
       a(k)=0.
      enddo
c  Loop over frames to compute mean frame
      do kb=1,np
       kp=(kb-1)*lb
       do k=1,lb
        ka=kp+k
        a(k)=a(k)+x(ka)
       enddo
      enddo
c------
      sm=0.d0
      sv=0.d0
      bk=float(np)
      do k=1,lb
       a(k)=a(k)/bk
       sm=sm+dble(a(k))
       sv=sv+dble(a(k)*a(k))
      enddo
      sm=sm/dfloat(lb)
      sg=sngl(dsqrt((sv-dfloat(lb)*sm*sm)/dfloat(lb-1)))
c  Subtract mean a from x
      do kb=1,np
       kp=(kb-1)*lb
       do k=1,lb
        ka=kp+k
        r(ka)=x(ka)-a(k)
       enddo
      enddo
c------
      return
      end
c
c===========================================
c
      subroutine spec(x,w,sp,f,v,n,lb,npd,il,jf,nt,sw,io)
c********************************************************
c  Frame averaged spectrum  M. J. Hinich 9-26-2005
c   x - data, w - taper weights
c   n - sample size, lb - frame size, nt - taper, nw - smooth sp
c   il - low freq bin of passband, iu - high freq bin of passband
c========================================================
      allocatable::y
      real x(n),sp(jf),w(lb),y(:)
      real*8 v(4*n+15)
      complex*16 f(0:n-1),cs2
      intent(in) x,w,n,lb,il,jf,nt,sw,io
      intent(out) sp
c**********************
      allocate(y(0:lb-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for array a'
       stop
      endif
c  Initialize arrays
      do k=1,jf
       sp(k)=0.
      enddo
c-------
      np=n/lb
      bk=float(np)
      rp=float(npd)
c  Scale factor
      cs2=dcmplx(dfloat(lb))
c  Pointers
      ks=il-1
      iu=ks+jf
c  Start loop over frames
      do kb=1,np
       kp=(kb-1)*lb
c  Taper frame kb
       do k=1,lb
        ka=kp+k
        if(nt>0) then
         y(k-1)=w(k)*x(ka)
          else
           y(k-1)=x(ka)
        endif
c  Move x frame to complex f
        f(k-1)=dcmplx(dble(y(k-1)))
       enddo
c  FFT of frame
       call cfftf(lb,f,v)
c-------
c  Periodogram fl<k<=fu
       do k1=il,iu
        jk=k1-ks
        sp(jk)=sp(jk)+sngl(dreal(f(k1)/cs2*f(lb-k1)))
        if(kb==np) then
         sp(jk)=sp(jk)/bk*rp
         if(nt>0) then
          sp(jk)=sp(jk)/sw
         endif
         if(sp(jk)<1.e-6) then
          sp(jk)=1.e-6
         endif
        endif
       enddo
      enddo
c  End of loop on frames
c-------
      deallocate(y)
      return
      end
c
c===========================================
c
      subroutine grams(a,gr,sp,w,th,lb,il,jf,seg,nsg,kp,io)
c********************************************
c  Input: a - chi**2, sp -spectrum, seg - segments, kp - segment count
c  Output: gr(1,:) - sigcoh probs, gr(jf+1:) - dB spectrum for column kr
c          w - array of signal coherences
c============================================
      real a(jf),sp(jf),w(jf),gr(2*jf,seg),ach,sch,th,thrs
      integer seg
      intent(in) a,sp,th,lb,il,jf,seg,nsg,kp,io
      intent(out) gr,w
c********************************************
      if(seg>1) then
c  Coherence threshold
       thrs=1.0-1.0/float(jf)
      endif
      if(seg==1) then
       thrs=th
      endif
      rl=float(lb)
      scc=float(nsg)/(rl*rl)
      ks=il-1
      do k=1,jf
       ach=a(k)/rl
       sch=sqrt(ach/(1.0+ach))
c  Signal coherence =>w
       if(sch>=thrs) then
c  Coherence is larger than threshold
        w(k)=sch
         else
          w(k)=0.0
       endif
       call cdchi((2.0*scc*a(k)),2.0,tpr,io)
c  Coherence probability
       gr(k,kp)=tpr
c  Spectrum in dB
       gr(jf+k,kp)=10.0*alog10(sp(k))
      enddo
c  End of frequencies loop
c-------
      return
      end
c
c===========================================
c
      subroutine cohsig(w,z,sg,f,g,v,lb,il,jf,sr,th,iunit,nc,seg,kr,io)
c********************************************
c  Input: w - signal coherences, v - initialized =>lb, th - coherence threshold
c  Output: z - CP of mean frame => column kr
c============================================
      real w(jf),z(lb,nc)
      real*8 v(4*lb+15),sm,sv
      complex*16 f(0:lb-1),g(0:lb-1)
      integer seg
      character*11 fry
      character*8 per
      character*2 chno,co
      intent(in) w,sg,g,v,lb,il,jf,sr,th,iunit,nc,seg,kr,io
      intent(out) z
c********************************************
c  Coherent average signal
      n2=lb/2
      ks=il-1
      iu=ks+jf
      kount=0
      fry='-frequency='
      per='-period='
      nc2=nc+nc+6
      rl=float(lb)
c  Resolution bandwidth
      if(iunit==1) then
       rb=sr/rl
      endif
      if(iunit==2) then
       rb=sr/rl*1.e3
      endif
c-------
      write(io,*)
c  Loop to zero outband & conjugate within band
      do k=1,n2
       if(k<il) then
        f(k)=(0.d0,0.d0)
        f(lb-k)=(0.d0,0.d0)
         elseif(iu<n2 .and. k>iu) then
          f(k)=(0.d0,0.d0)
          f(lb-k)=(0.d0,0.d0)
           else
            tpr=w(k-ks)
c  Test coherence against threhold
            if(tpr>th) then
             f(k)=tpr*dconjg(g(k))
             f(lb-k)=dconjg(f(k))
             kount=kount+1
             if(seg==1) then
c  Output frequency or period
              if(iunit>0) then
               fk=float(k)*rb
               write(io,'(a2,a11,f12.5)')chno(kount),fry,fk
              endif
              if(iunit==0) then
               fk=(rl*sr)/float(k)
               write(io,'(a2,a8,f12.5)')chno(kount),per,fk
              endif
             endif
              else
               f(k)=(0.d0,0.d0)
               f(lb-k)=(0.d0,0.d0)
            endif
       endif
      enddo
c  Trap no coherences >th
      if(kount==0) then
       return
      endif
c  Mean
      f(0)=g(0)
      call cfftf(lb,f,v)
c  Rescale signal
      sm=0.d0
      sv=0.d0
      do k=1,lb
       st=real(f(k-1))/float(lb)
       z(k,kr)=st
       sm=sm+dble(st)
       sv=sv+dble(st*st)
      enddo
      sm=sm/dfloat(lb)
      sig=sngl(dsqrt((sv-dfloat(lb)*sm*sm)/dfloat(lb-1)))
c  Scaling ratio
      scal=sg/sig
      do k=1,lb
       st=scal*z(k,kr)
       z(k,kr)=st
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine datc(x,ts,fmt,mk,nmk,is,ie,n,nc,ii,io)
c********************************************************
c  Reads n data values from nc columns
c   x - input data matrix, ts - time marks if indicated by 'm'
c========================================================
      real x(n,nc)
      character*50 fmt
      character*25 ts(nmk)
      character mk
      intent(in) fmt,mk,nmk,is,ie,n,nc,ii,io
      intent(out) x,ts
c********************************************
      nt=0
      if(mk=='m') then
       ia=index(fmt,'(')
       ib=index(fmt,'a')
       if(ib==ia+1) then
c  Time mark column on left
        nt=1
         elseif(ib>ia+2) then
c  Time mark column on right
          nt=2
           elseif(nt==0) then
            write(io,'(7x,''/Data & time format is incorrect''/)')
            write(io,'(a50)')fmt
            write(6,'(7x,''/Data & time format is incorrect''/)')
            write(6,'(a50)')fmt
            write(109,'(7x,''/Data & time format is incorrect''/)')
            write(109,'(a50)')fmt
            stop
       endif
      endif
c-------
      rewind(ii)
c  Skip header and parameter records
      read(ii,*)
      read(ii,*)
c  Read all the data
      if(mk=='m') then
       if(is==1) then
        do i=1,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)ts(i),(x(i,k),k=1,nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc),ts(i)
         endif
        enddo
       endif
       if(is>1) then
c  Read is:ie
        do i=1,is-1
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         if(nt==1) then
          read(ii,fmt,end=1,err=2,iostat=ko)ts(i),(x(i-is+1,k),k=1,nc)
         endif
         if(nt==2) then
          read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc),ts(i)
         endif
        enddo
       endif
      endif
      if(mk/='m') then
       if(is==1) then
        do i=1,n
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i,k),k=1,nc)
        enddo
       endif
       if(is>1) then
        do i=1,is-1
         read(ii,fmt,end=1,err=2,iostat=ko)
        enddo
        do i=is,ie
         read(ii,fmt,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc)
        enddo
       endif
      endif
      return
c--------
  1   write(io,*)'**** eof in read from data file '
      write(io,*)'i = ',i
      return
      write(6,*)'**** eof in read from data file '
      write(6,*)'i = ',i
      write(109,*)'**** eof in read from data file '
      write(109,*)'i = ',i
c  Error in data read
  2   write(io,*)'** Error in read of data file ',ko
      write(io,*)'Check format',fmt
      write(6,*)'** Error in read of data file ',ko
      write(6,*)'Check format',fmt
      write(109,*)'** Error in read of data file ',ko
      write(109,*)'Check format',fmt
      stop
      end
c
c===========================================
c
      subroutine interp(y,w,f,v,n,npd,io)
c********************************************************
c  Interpolate by padding FFT with 0's and then inverse FFT for data
c========================================================
      real y(npd*n),w(npd*n)
      real*8 v(4*npd*n+15)
      integer t
      complex*16 f(0:npd*n-1)
      intent(in) y,n,npd,io
      intent(out) w
c**************************
      if(npd<2) then
       return
      endif
c--------
      no=npd*n
      n2=n/2
      no2=no/2
c  pi/2
      pi2=acos(-1.)/2.
c  No. tapered after n2
      ntap=max(nint(.05*float(n)),10)
c  Move x to dcomplex f
      do t=1,n
       f(t-1)=dcmplx(dble(y(t)))
      enddo
c-------
c  Initialize v
      call cffti(n,v)
c  FFT of x
      call cfftf(n,f,v)
      do k=1,n2
       f(no-k)=dconjg(f(k))
      enddo
      do k=n2+1,n2+ntap
       f(k)=cos(pi2/float(ntap+1)*float(k-n2))*f(k)
       f(no-k)=dconjg(f(k))
      enddo
      do k=n2+ntap,no2
       f(k)=(0.d0,0.d0)
       f(no-k)=(0.d0,0.d0)
      enddo
      call cffti(no,v)
      call cfftf(no,f,v)
      do t=1,no
       w(t)=real(f(t-1))/float(no)
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine wrh(su,sr,rb,n,ov,nsg,seg,npd,lb,pt,p,pv,th,nw,fl,fu,iu
     *nit,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      real*8 dlb
      integer p,ov,seg
      character*10 su
      intent(in) su,sr,rb,ov,nsg,seg,npd,lb,pt,p,pv,th,nw,fl,fu,iunit,io
c**********************
      dlb=dfloat(lb)
      il=nint(sngl(dble(fl)*dlb))
      il=max(il,1)
      iu=nint(sngl(dble(fu)*dlb))
      iu=min(iu,n/2)
      rl=float(lb)
      rpd=float(npd)
      write(io,'(60(''='')/)')
      write(io,'('' Overlap = '',i5,'' Segment Length = '',i7,'' Number
     *of Segments = '',i6/)')ov,nsg,seg
      if(iunit==1) then
       a=float(il)/rl*sr
       b=float(iu)/rl*sr
c  Lower and upper frequencies for band
       write(io,'('' Sampling rate ='',f12.3,'' kHz''/)')sr/rpd
       write(io,'('' Resolution Bandwidth = '',f10.3,'' Hz.'',2x,''No. F
     *rames = '',i7,1x,''Frame Length = '',i7/)')rb,n/lb,lb/npd
       write(io,'('' Passband ('',f9.3,2x,f10.3,'' ) kHz'')')a,b
       if(npd>1) then
        write(io,'(/7x,''Data interpolated by a factor of '',i5)')npd
       endif
      endif
c------
      if(iunit==2) then
       a=float(il)/rl*sr*1.e3
       b=float(iu)/rl*sr*1.e3
       write(io,'('' Sampling rate ='',f12.3,'' Hz''/)')1.e3*sr/rpd
       write(io,'('' Resolution Bandwidth = '',f10.5,'' Hz.'',2x,''No. F
     *rames = '',i7,1x,''Frame Length = '',i7/)')rb,n/lb,lb/npd
       write(io,'('' Passband ('',f7.2,2x,f9.2,'' ) Hz'')')a,b
       if(npd>1) then
        write(io,'(/7x,''Data interpolated by a factor of '',i5)')npd
       endif
      endif
c------
      if(iunit==0) then
       a=rb/float(il)
       b=rb/float(iu)
       write(io,'('' Sampling interval = '',g12.6,1x,a7/)')sr*rpd,su
       write(io,'('' Resolution Bandwidth = '',g10.4,1x,a7,2x,''No. Fram
     *es = '',i7,1x,''Frame Length = '',i7/)')rb,su,n/lb,lb/npd
       write(io,'('' Passband ('',f9.2,2x,f9.2,'' )'',1x,a7)')a,b,su
       if(npd>1) then
        write(io,'(/7x,''Data interpolated by a factor of '',i3)')npd
       endif
      endif
c-------
      write(io,'(/7x,''No. of frequencies in band '',i6/)')iu-il+1
c  Number of values tapered from each end
      nt=nint(pt*rl/200.)
      if(pt>0.0.and.nt>=1) then
       write(io,'('' % taper = '',f5.2,4x,''No. tapered at each end = ''
     *,i5/)') pt,nt
      endif
      if(pt>0.0.and.nt<1) then
       write(io,'('' % Taper yields no. of tapered values < 1''/)')
      endif
c-------
      if(nw>=3) then
       kw=2*(nw/2)+1
       write(io,'(2x,''Spectrum is smoothed by a truncated cosine using
     *'',i5,'' adjacent bins''/)') kw
        elseif(nw>0.and.nw<3) then
         write(io,'('' Smoothing bandwidth < resolution bandwidth'')')
         write(io,'('' Spectrum is not smoothed''/)')
          else
          write(io,'('' Spectrum is not smoothed''/)')
      endif
      write(io,'(2x,''Coherence threshold for the coherent part of the f
     *rame = '',f7.4/)')th
      if(p>0) then
       write(io,'(2x,''Probability tail value threshold for the AR fit i
     *terations = '',f7.4/)')pv
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine wrstat(i,am,sig,sk,c4,c6,smax,smin)
c****************************************************
c  Data statistics
c======================
      implicit none
      real am,sig,sk,c4,c6,smax,smin
      integer i
      intent(in) i,am,sig,sk,c4,c6,smax,smin
c********************************************
      write(i,'(60(''='')/)')
      write(i,'('' Mean = '',g12.3,7x,''Std Dev = '',g12.3/)')am,sig
      write(i,'('' Skew = '',g10.3,4x,''Kurtosis = '',g10.3,4x,''C(6) =
     *'',f16.3/)')sk,c4,c6
      write(i,'('' Max value = '',g12.3,7x,''Min value = '',g12.3)')smax
     *,smin
c-------
      return
      end
c
c===========================================
c
      subroutine wrtspec(gr,w,b,lb,il,jf,seg,su,sr,iunit,m,io)
c*******************************************
c  Output to file m - ,coherence w, coprob gr(1:), spectrum gr(fk+1:)
c  b - periodogram of mean frame
c===========================================
      real gr(2*jf,seg),w(jf),b(jf)
      integer seg
      character*10 su
      intent(in) gr,b,w,lb,il,jf,seg,su,sr,iunit,m,io
c*************************
      rl=float(lb)
      ks=il-1
      if(iunit==1) then
c  Resolution bandwidth
       rb=sr/rl
       do k=1,jf
        fk=float(ks+k)*rb
        spc=gr(jf+k,1)
        write(m,'(f14.4,2x,f10.4,2x,f5.3,2x,f7.5,2x,f16.3)')fk,spc,w(k),
     *gr(k,1),b(k)
       enddo
      endif
      if(iunit==2) then
       rb=sr/rl*1.e3
       do k=1,jf
        fk=float(ks+k)*rb
        spc=gr(jf+k,1)
        write(m,'(f14.4,2x,f10.4,2x,f5.3,2x,f7.5,2x,e16.4)')fk,spc,w(k),
     *gr(k,1),b(k)
       enddo
      endif
      if(iunit==0) then
c  Periods
       do k=1,jf
        rk=(rl*sr)/float(ks+k)
        spc=gr(jf+k,1)
        write(m,'(f14.4,2x,f10.4,2x,f5.3,2x,f7.5,2x,e16.4)')rk,spc,w(k),
     *gr(k,1),b(k)
       enddo
      endif
      call flush(m)
c-------
      return
      end
c
c===========================================
c
      subroutine wrtgrmc(x,sr,su,lb,il,jf,seg,iunit,nc,kr,km,io)
c********************************************
c  Write spectrogram & sigcohgram for each segment kr
c  x(1,seg) - sigcoh probabilities, x(jf+1,seg) - dB spectrum (gr in call)
c  sr - sampling rate, jf - no. spectral values, seg - segments, nc - columns
c  kr - columnindex, km - segment start with or without overlap
c  gr(1,:) - sigcoh probs, gr(jf+1:) - spectrum for column kr
c  M. J. Hinich   Version: 11-14-2009
c===========================================================
      allocatable::frq
      real x(2*jf,seg),frq(:),start
      character*10 su
      integer c,seg,lb,il,jf,iunit,nc,kr,mk,io
      intent(in) x,sr,su,lb,il,jf,seg,iunit,nc,kr,km,io
c********************************************
c  Segment start time in minutes
      if(iunit==1) then
       start=float(km)/(60.0*sr)
        elseif(iunit==2) then
         start=float(km)/(60.0*sr*1.e3)
          elseif(iunit==0) then
           start=float(km)
      endif
c-------
      n1=nc+6+kr
      n2=nc+nc+6+kr
      allocate(frq(0:jf-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for array a'
       stop
      endif
      rl=float(lb)
      do k=0,jf-1
c  Frequencies or periods
       if(iunit==1) then
        frq(k)=float(il+k)*(sr/rl)
         elseif(iunit==2) then
          frq(k)=float(il+k)*(sr/rl*1.e3)
           elseif(iunit==0) then
            frq(k)=(rl*sr)/float(il+k)
       endif
      enddo
c  Spectrogram
      do c=1,seg
       write(n1,'(:,10000(f14.5,1x))')(x(jf+k,c),k=1,jf)
c  Sigcohgram
       write(n2,'(:,10000(f14.5,1x))')(x(k,c),k=1,jf)
      enddo
      write(n1,'(:,10000(f14.5,1x))')(frq(k),k=0,jf-1)
      write(n2,'(:,10000(f14.5,1x))')(frq(k),k=0,jf-1)
c-------
      deallocate(frq)
      return
      end
c
c===========================================
c
      subroutine rmaxmin(x,jf,maxx,minx,amin,amax)
c********************************************
c  Maximum & minimum values of array x
c===========================================================
      implicit none
      real x(jf),maxx,minx,amin,amax
      integer jf,t
      intent(in) x,jf,amin,amax
      intent(out) maxx,minx
c********************************************
         maxx=amax
         minx=amin
         do t=1,jf
          if(x(t)>maxx) then
           maxx=x(t)
          endif
          if(x(t)<minx) then
           minx=x(t)
          endif
         enddo
c-------
      return
      end
c===========================================
c
      function chno(kount)
c********************************************
c  Convert integer kount<100 to charactger*2 chno
c  M. J. Hinich   Version: 7-15-2007
c===========================================================
      integer kount
      character*2 chno
      intent(in) kount
c********************************************
       if(kount<10) then
        chno(1:1)=achar(kount+48)
        chno(2:2)=' '
       endif
       if(kount>=10.and.kount<20) then
        chno(1:1)='1'
        chno(2:2)=achar(kount-10+48)
       endif
       if(kount>=20.and.kount<30) then
        chno(1:1)='2'
        chno(2:2)=achar(kount-20+48)
       endif
       if(kount>=30.and.kount<40) then
        chno(1:1)='3'
        chno(2:2)=achar(kount-30+48)
       endif
       if(kount>=40.and.kount<50) then
        chno(1:1)='4'
        chno(2:2)=achar(kount-40+48)
       endif
       if(kount>=50.and.kount<60) then
        chno(1:1)='5'
        chno(2:2)=achar(kount-50+48)
       endif
       if(kount>=60.and.kount<70) then
        chno(1:1)='6'
        chno(2:2)=achar(kount-60+48)
       endif
       if(kount>=70.and.kount<80) then
        chno(1:1)='7'
        chno(2:2)=achar(kount-70+48)
       endif
       if(kount>=80.and.kount<90) then
        chno(1:1)='8'
        chno(2:2)=achar(kount-80+48)
       endif
       if(kount>=90.and.kount<100) then
        chno(1:1)='9'
        chno(2:2)=achar(kount-90+48)
       endif
c-------
      return
      end