program trispec
c*********************************************************************
c  Version 8-17-2008   Melvin Hinich
c  File 3 (fname.out) is always opened where 'fname' is the name of the data
c  file used to the left of the .extention. The output in 3 are statistics for
c  the trispectrum for the data used along with summary stats for this data.
c  IMPORTANT - Error messages are written to trispec-error.out.
c
c   File 2 - data file, 7 - k3 > 0 trispectrum file, 4 - k3 < 0 trispectrum
c  The trispectrum file name is #.trispectrum where # = data file name.
c======================================================================
c  Parameters in control file RUNTRISPEC.CNL. Use spaces between entries!
c  A delimiter symbol * must be placed on each line.
c  Place comments on the lines of runtrispec.cnl file after the symbol *.
c  Runtrispec.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  Lines for the control file RUNTRISPEC.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 runtrispec.cnl. The integer to the right
c  of the = sign in runtrispec.cnl is read by the program and sets the number
c  of runs.
c  In any *.cnl called by runtrispec.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  IMPORTANT: Leave at least ONE SPACE after the paramter
c
c  ALWAYS use 1.0 rather than .0 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.
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  The number of runs is followed by a list of the path & names of the cnl
c  files that pass parameters to the program
c  For example if !nrun=2:
c  \prog\trispec-fil1.cnl             * cnl file
c  \prog\trispec-fil2.cnl             * cnl file
c
c  There has to be !nrun cnl files in runtrispec.cnl
c
c  These entries are followed by the line
c  # * End of Runtrispec.cnl control file names read
c
c  The program seeks the delimiter symbol # to end the reading of the
c  control file names.
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 control file reads
c
c  The program seaks the delimiter symbol # to end the reading of the
c  parameters for each control file set in Runtrispec.cnl
c
c=========================
c
c  Data file name=data file name with path
c
c=========================
c
c  Use column number=!ncol
c  !ncol is the column number of the nc columns in the data file
c  Set !ncol=1 if there is only one column in the data (flat file)
c
c=========================
c
c  Resolution bandwidth=!rb  ( !r - real wildcard)
c   This bandwidth is in Hz or the unit used below if frequency is not typed
c  Note: rb determines the frame length lb used to compute the spectrum
c
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  NOTE: sr, the 2nd number on line 2 of the data file is the sampling
c   rate in kHz if frequency+K above, or Hz if frequency+H above, or sr is
c   a sampling interval with the unit determined by the unit determined by
c   the characters placed after unit=
c   If there is an 'a' in the format, then a time mark file is read
c   An example string is: (f7.2,2x,a7)
c   The time mark must be less than 21 characters long
c
c=======================================================================
c
c  Lower frequency=!rl  upper frequency=!ru or
c  Longest period=!rl   shortest period=!ru
c   The passband is (!rl,!ru)
c   !rl - lower frequency, !ru - upper frequency of band in units of kHz
c   if frequency+kHz  is typed or if OR Hz if frequency+Hz is typed on a line
c   !rl - longest period - the longest period
c   !ru - shortest period - the shortest period if a unit is used
c
c   If !ru= 0.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   !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  % taper=!rpt
c   !rpt = % taper of frame for sidelobe reduction (0.0 < rpt < 25.)
c   If !rpt = 0.0, frames are NOT tapered
c
c=========================
c
c  Number of resamples=!n, quantile=!q
c   Use a positive integer !n > 10 for the number of resamples if you want
c   to bootstrap the nonlinearity test statistic for the !q'th
c   quantile of the normalized trispectrum.
c   The bootstrap used is a uniform shuffle of the residuals of the
c   optimizing AR(p) fit used to whiten the data. The starting value of
c   p is set in %.cnl.
c   Type 0 if you do not want to bootstrap the test statistics
c   Note - bootstrapped p - value fractiles are only computed when the data
c   is whitened with an AR fit.
c
c=========================
c
c  NOTE:
c  The modified Hinich (1982) nonlinearity test uses the quantile qv = 0.8
c  Now I prefer qv = 0.99
c
c=========================
c
c  If you want to prewhiten with an AR(np) fit type the character string
c   prewhiten with AR(!np) fit
c   !np - starting number of parameters of a recursive least squares AR
c    fit if np > 0.0 If np = 0 then this step is bypassed.
c  No AR fit is made if this character string is absent in the lines above #.
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   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  Prob threshold=!rpv (for AR t values if AR fit is used)
c   !rpv - probability level threshold for the t statistics for the AR
c    coefficients of each AR fit. If a coefficient has a t value whose
c    absolute value has a probability Pr(abs(t)) < 1 - rpv, it is deleted in
c    the next AR fit. The covariance matrix is adjusted for the subset of
c    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  Data file format for the first two lines (records)
c
c  Line 1: A character string that identifies the data in thie file .
c  Line 2: If the unit is in frequency then line 2 using an EXAMPLE format is
c  rows=!nfl, columns=!nc, sampling rate=!sr, format=(f10.5,2x,a12) unit=kHz
c  Line 2: If the unit is a time unit then line 2 is
c  rows=!nfl, columns=!nc, sampling interval=!sr, format=(f10.5,2x,a12) unit=day
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.
c  Up to 50 columns can be read.
c
c   nfl - no. of observations, sr - sampling rate in kHz if freq. is
c   set, or in multiples of the time unit set in trispec-fn.cnl
c
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   This character string is used to check the read from %.cnl
c
c   If there is an 'a' in the format then a time & data character string
c   is read. The time & data string must be less than 21 characters long.
c   An example format for a data file with two columns and a date & time
c   column 7 spaces wide is: (2(f7.2,1x),2x,a7).
c
c*****************************************************
      integer run,nr
      real gs
      character*700 name
      character*80 buf,cnl
      character*50 charead
      character*20 par
      character*10 su
      character*2 delim
      logical exs
c**********************
c  Start time
      call cpu_time(gs)
      open(9,file='Trispec-error.out')
c  Inquire if file exists
      inquire(file='Runtrispec.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'Runtrispec.cnl file does not exist'
       write(9,*)'Runstrispec.cnl file does not exist'
       stop
      endif
c  Open Runtrispec.cnl
      open(1,file='Runtrispec.cnl',err=1,iostat=io1,status='old')
      read(1,'(a80)',end=2,err=3,iostat=io2) buf
      ig=index(buf,'*')
      if(ig==0) then
       write(6,*)'Place comment marker * after each line in Runbispect.c
     *nl'
       write(9,*)'Place comment marker * after each line in Runbispect.c
     *nl'
       write(6,'(a70)')buf
       write(9,'(a70)')buf
       stop
      endif
c  No. of runs
      par='uns'
      delim(1:2)='= '
      run=numread(par,buf,delim,ia,ier,io)
c  Multiple runs
      do nr=1,run
c  Read control file
       read(1,'(a80)',end=4,err=5,iostat=io2) cnl
       ig=index(cnl,'*')
       cnl=cnl(:ig-1)
       inquire(file=cnl,exist=exs)
       if(.not.exs) then
        write(6,'(/)')
        write(6,*) 'cnl file does not exist ',cnl
        write(9,'(/)')
        write(9,*) 'cnl file does not exist ',cnl
        stop
       endif
c  Open control file
       open(11,file=cnl,err=6,iostat=io3,status='old')
       call top(gs,run,nr)
      enddo
c  End of loop on runs
      stop
  1   write(6,*) '**** Error on opening Runtrispect.cnl ',io1
      stop
  2   write(6,*) 'End of file on runtrispec.cnl reads'
      write(9,*) 'End of file on runtrispec.cnl reads'
  3   write(6,*) 'Place the symbol # after the two top reads ',io1
      write(9,*) 'Place the symbol # after the two top reads ',io1
      stop
  4   write(6,*) 'End of file on runtrispec.cnl reads'
      write(9,*) 'End of file on runtripsec.cnl reads'
      stop
  5   write(6,*) '**** Error on reading name of cnl file ',nr,io2
      write(6,*) cnl
      write(9,*) '**** Error on reading name of cnl file ',nr,io2
      write(9,*) cnl
      stop
  6   write(6,*) '**** Error on opening cnl file no. ',nr,io3
      write(6,*) cnl
      write(9,*) '**** Error on opening cnl file no. ',nr,io3
      write(9,*) cnl
      stop
c-------
      stop
      end
c
c===========================================
c
      subroutine top(gs,run,nr)
c********************************************************
c  Open data file and read data
c========================================================
      parameter(mt=10,ndy=26,ny=2009)
      allocatable::a,b,c,f,in,r,rt,s,sp,tr,ti,t1,t2,t3,u,v,x,w,z
      real*8 c(:),sp(:),tr(:),ti(:),v(:),b(:),x(:,:),r(:),s(:),u(:,:,:),
     *w(:),z(:,:,:),rlb
      real qv
      complex a(:),rt(:)
      complex*16 f(:)
      integer t1(:),t2(:),t3(:),in(:),ot,p,po,q,res,run
      character*700 name
      character*80 buf,id
      character*50 charead,fmt,fname
      character*24 dt
      character*20 par
      character*15 su
      character*4 str
      character*2 delim
      logical exs
      intent(in) gs,run,nr
c********************************************
      ip=0
      do k=1,17
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(11,'(a80)') buf
       ia=index(buf,'#')
       if(ia>0) then
        exit
       endif
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after parameters on map.cnl'
        write(6,'(a70)')buf
        write(9,*)'Place comment marker * after parameters on map.cnl'
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
c  Data file name from control file
      par='ile name'
      delim(1:2)='= '
      fmt=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(6,*)'Error in reading the data file name from %.cnl'
       write(9,*)'Error in reading the data file name from %.cnl'
       write(6,*)'Error number for charead = ',ier
       write(9,*)'Error number for charead = ',ier
       stop
      endif
c-------
      close(2)
c  Inquire if file exists
      inquire(file=fmt,exist=exs)
      if(.not.exs) then
       write(6,*)'Input file ',fmt,' does not exist'
       write(9,*)'Input file ',fmt,' does not exist'
       stop
      endif
      open(2,file=fmt,err=1,iostat=io1,status='old')
c-------
c  Trispectrum file
      ia=index(fmt,'.')
      if(ia==0) then
       write(6,*)'Data file name does not have an extension'
       write(6,*)fmt
       write(9,*)'Data file name does not have an extension'
       write(9,*)fmt
       stop
      endif
c  Parse for up to 5 path switches
      ib=index(fmt,'\')
      if(ib>0) then
       ic=index(fmt(ib:),'\')
       ib=ib+ic
       if(ic>0) then
        ic=index(fmt(ib:),'\')
        ib=ib+ic
        if(ic>0) then
         ic=index(fmt(ib:),'\')
         ib=ib+ic
        endif
        if(ic>0) then
         ic=index(fmt(ib:),'\')
         ib=ib+ic
        endif
       endif
        else
         ib=1
      endif
      fname=fmt(ib:ia)
      nfn=ia-ib+2
      fname(nfn-1:nfn+11)='_Trispec.out'
c  Open summary output file
      close(3)
      open(3,file=fname,err=4,iostat=io4)
c  Write name of test
      write(3,'(11x,'' Trispectrum Statistics by M. J. Hinich'')')
      call clock(gs,0,3)
      fname(nfn-1:nfn+15)='-k3+_trispec.out'
c  Open trispectrum output files for graphing
      close(7)
      open(7,file=fname,err=2,iostat=io2)
      call fdate(dt)
      write(7,'('' Version '',i2,''-'',i2,''-'',i4,7x,''Run '',a24)') mt
     *,ndy,ny,dt
      fname(nfn-1:nfn+15)='-k3-_trispec.out'
      close(4)
      open(4,file=fname,err=3,iostat=io3)
      call fdate(dt)
      write(4,'('' Version '',i2,''-'',i2,''-'',i4,7x,''Run '',a24)')mt,
     *ndy,ny,dt
c-------
c  Column number
       kcol=0
       par='umber'
       kcol=numread(par,name,delim,ia,ier,io)
       if(ier>0.or.kcol<=0) then
        write(3,'(7x,''ERROR reading use column number '',i3)')kcol
        write(6,'(7x,''ERROR reading use column number '',i3)')kcol
        write(9,'(7x,''ERROR reading use column number '',i3)')kcol
        stop
       endif
c-------
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(9,'(/'' 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
       iunit=0
      endif
c  Trap iunit
      if(iunit==7) then
       write(6,'(/'' Sampling unit not set in the cnl file'',a15)')buf
       write(9,'(/'' Sampling unit not set in the cnl file'',a15)')buf
       stop
      endif
c-------
      rewind(2)
c  Read data file header
      call readhead(id,nfl,nc,sr,fmt,iunit,2,3)
      write(3,'(/7x,''Column Number Used = '',i3,2x,''No. of Columns = '
     *',i3)')kcol,nc
      if(iunit>0) then
c  Lower frequency
       par='wer freq'
       fl=rnumread(par,name,delim,ia,ier,io)
       if(ier==1) then
        write(3,'(/'' Type lower frequency rather than longest period f
     *or the passband in %.cnl'')')
        write(6,'(/'' Type lower frequency rather than longest period fo
     *r the passband in %.cnl'')')
        stop
       endif
c  Upper frequency
       par='per freq'
       fu=rnumread(par,name,delim,ia,ier,io)
       if(ier==1) then
        write(3,'(/'' Type upper frequency rather than shortest period f
     *or the passband in %.cnl'')')
        write(6,'(/'' Type upper frequency rather than shortest period f
     *or the passband in %.cnl'')')
        stop
       endif
        else
c  Longest period
         par='gest period'
         fl=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(3,'(/'' Type longest period rather than lower frequency
     *for the passband in %.cnl'')')
          write(6,'(/'' Type longest period rather than lower frequency
     *for the passband in %.cnl'')')
          stop
         endif
c  Shortest period
         par='test period'
         fu=rnumread(par,name,delim,ia,ier,io)
         if(ier==1) then
          write(3,'(/'' Type shortest period rather than upper frequency
     * for the passband in %.cnl'')')
          write(6,'(/'' Type shortest period rather than upper frequency
     * for the passband in %.cnl'')')
          stop
         endif
      endif
c  % taper
      par='taper'
      pt=rnumread(par,name,delim,ia,ier,io)
c  Resolution bandwidth
      par='idth'
      rb=rnumread(par,name,delim,ia,ier,io)
      if(ier==4) then
       write(3,*)'Place a period in the resolution bandwidth value'
       ia=index(name,'idth')
       write(3,*)name(ia+5:ia+15)
       stop
      endif
      if(rb<=0.0) then
       write(3,*)'Resolution bandwidth < 0'
       stop
      endif
      p=0
c  Prewhiten by an AR(p) fit
      ia=index(name,'rewhite')
      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,delim,ia,ier,io)
         elseif(ic>0) then
          par='r'
          delim(1:2)='()'
          p=numread(par,name,delim,ia,ier,io)
       endif
      endif
      delim(1:2)='= '
      if(p>0) then
c  Prob threshold for AR t values
       par='hreshold'
       pv=rnumread(par,name,delim,ia,ier,io)
c  Trap pv
       if(pv<=0.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'
        stop
       endif
       if(pv>1) then
        write(3,*)'Prob level pv ',pv,' > 1'
        write(6,*)'Prob level pv ',pv,' > 1'
        stop
       endif
      endif
c-------
      if(iunit==0) then
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
c-------
      if(fu>0.0 .and. fl>fu) then
       write(3,*)'Upper frequency of band ',fu,' < lower frequency ',fl
       write(6,*)'Upper frequency of band ',fu,' < lower frequency ',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'
         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'
           else
            write(3,*)'Check to see the (fl,fu) units are ',su
            write(6,*)'Check to see the (fl,fu) units are ',su
       endif
       stop
      endif
c-------
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
       stop
      endif
c  File read start
      par='start'
      is=numread(par,name,delim,ia,ier,io)
      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
       stop
      endif
c  File read end
      par='end'
      ie=numread(par,name,delim,ia,ier,io)
      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
       stop
      endif
      if(is<1) then
       write(3,*)'IS in CNL file is < 1 ',is
       write(6,*)'IS in CNL file is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(3,*)'IE < IS in CNL file ',is,ie
         write(6,*)'IE < IS in CNL file ',is,ie
         stop
      endif
c-------
c  Resamples
      res=0
      par='sample'
      res=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of resamples'
       write(3,*)par
       write(6,*)'Error ',ier,' in read of resamples'
       write(6,*)par
       write(9,*)'Error ',ier,' in read of resamples'
       write(9,*)par
       stop
      endif
c  Trap res
      if(res>0.and.res<11) then
       write(3,'(/2x,''Bootstrap is turned off if the number of resampl
     *s is less than 11''/,'' Resamples = '',i2)')res
       write(6,'(/2x,''Bootstrap is turned off if the number of resample
     *s is less than 11''/,'' Resamples = '',i2)') res
       write(9,'(/2x,''Bootstrap is turned off if the number of resample
     *s is less than 11''/,'' Resamples = '',i2)') res
       res=0
      endif
      if(p==0) then
       res=0
      endif
c  Quantile
      qv=0.9
      par='ntile'
      qv=rnumread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,*)'Error ',ier,' in read of quantile'
       write(6,*)'Error ',ier,' in read of quantile'
       write(9,*)'Error ',ier,' in read of quantile'
       stop
      endif
c--------
c  Trap top band
      if(fu>sr/2.) then
       write(3,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       write(6,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       write(9,*)'Upper frequency = ',fu,' > Nyquist freq. ',sr/2.
       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(9,*)'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(9,*)'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(9,*)'Check to see the (fl,fu) units are ',su
       endif
       write(3,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(3,*)'iunit = ',iunit
       write(6,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(6,*)'iunit = ',iunit
       write(9,*)' Sampling unit ',su,' Sampling rate/interval ',sr
       write(9,*)'iunit = ',iunit
       stop
      endif
c  Convert to normalized values (0,.5)
      if(iunit>0) then
       gl=fl/sr
       if(fu>0.0) then
        gu=fu/sr
         else
          gu=0.5
       endif
        else
         gl=fl*sr
         gu=fu*sr
      endif
c-------
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
       stop
        elseif(ie==0) then
         ie=nfl
      endif
      n=ie-is+1
c  Find frame length lb
      if(iunit==1) then
c  kHz
       lb=nint(sr*1.e3/rb)
        elseif(iunit==2) then
c  Hz
         lb=nint(sr/rb)
          else
           lb=nint(rb/sr)
      endif
c-------
      if(lb<=9) then
       write(3,*)'Frame length = ',lb,' < 10'
       write(3,*)'Decrease resolution bandwidth = ',rb
       write(6,*)'Frame length = ',lb,' < 10'
       write(6,*)'Decrease resolution bandwidth = ',rb
       write(9,*)'Frame length = ',lb,' < 10'
       write(9,*)'Decrease resolution bandwidth = ',rb
       stop
c  Trap n<lb
        elseif(n/lb<4) then
         write(3,*)'Frame length > n/4'
         write(3,*)'Increase resolution bandwidth = ',rb
         write(6,*)'Frame length > n/4'
         write(6,*)'Increase resolution bandwidth = ',rb
         write(9,*)'Frame length > n/4'
         write(9,*)'Increase resolution bandwidth = ',rb
         stop
      endif
c-------
c  Number of values tapered from each end
      nt=nint(pt*float(lb)/200.)
c  If no. of tapered values at ends < 2 set pt=0
      if(nt<2) then
       pt=0.0
      endif
      rlb=dfloat(lb)
      n2=lb/2
c  Set il - lower band integer limit
      il=nint(sngl(dble(gl)*rlb))
      il=max(il,1)
c  Set iu - upper band integer limit
      if(fu>0.0) then
       iu=nint(sngl(dble(gu)*rlb))
       iu=min(iu,n2)
        else
         iu=n2
      endif
      if(il==iu) then
       write(3,*)'fu-fl is too small. Increase fu-fl'
       write(6,*)'fu-fl is too small. Increase fu-fl'
       write(9,*)'fu-fl is too small. Increase fu-fl'
       stop
        elseif(il==iu .and. iu==n2) then
         write(3,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
         write(6,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
         write(9,*)'fl ',fl,'=',fu,' fu: ERROR in band settings in cnl'
        stop
      endif
c-------
c  No. of frequencies used
      jf=iu-il+1
      e=alog(float(lb))/alog(float(n))
      if(e<=0.01.or.e>=.99) then
       write(3,*)'e = ',e,' is out of bound'
       write(6,*)'e = ',e,' is out of bound'
       write(9,*)'e = ',e,' is out of bound'
       stop
        elseif(e>0.3333) then
         write(3,'(/2x,''Condition for Consistency of Estimates is Viola
     *ted - e ='',f7.5)')e
         write(3,'(7x,''Resolution Bandwidth = '',f10.2)')rb
         write(6,'(/2x,''Condition for Consistency of Estimates is Viola
     *ted - e ='',f7.5)')e
         write(6,'(7x,''Resolution Bandwidth = '',f10.2)')rb
      endif
c-------
c  Size arrays, npt - pd for k3>0, nnt -pd for k3<0, ot -aliasing set
      call size(il,iu,lb,npt,nnt,ot,k3u,k3d,io)
c  Trap number of trifrequncies in IT region
      nit=npt+nnt
      if(nit<=5) then
       write(3,*)'Number of trifrequencies to be computed ',nit,' < 5'
       write(3,*)'Increase the passband'
       write(6,*)'Number of trifrequencies to be computed ',nit,' < 5'
       write(6,*)'Increase the passband'
       write(9,*)'Number of trifrequencies to be computed ',nit,' < 5'
       write(9,*)'Increase the passband'
       stop
      endif
      write(3,'(/3x,''Stationary Trispectrum Values for k3 > 0 = '',i9)'
     *)npt
      if(npt==0) then
       write(3,'(/7x,''ERROR''/)')
       write(3,'(3x,''Stationary Trispectrum Values for k3 > 0 = '',i1)'
     *)npt
       write(3,'(3x,''Increase exponent exp to increase frame'',i5)')lb
       write(6,'(/7x,''ERROR''/)')
       write(6,'(3x,''Stationary Trispectrum Values for k3 > 0 = '',i1)'
     *)npt
       write(6,'(3x,''Increase exponent exp to increase frame'',i5)')lb
       write(9,'(/7x,''ERROR''/)')
       write(9,'(3x,''Stationary Trispectrum Values for k3 > 0 = '',i1)'
     *)npt
       write(9,'(3x,''Increase exponent exp to increase frame'',i5)')lb
       stop
      endif
      write(3,'(3x,''Stationary Trispectrum Values for k3 < 0 = '',i9)'
     *)nnt
      write(3,'(3x,''Total  Number of Values = '',i9)')npt+nnt
      write(3,'(3x,''Aliasing Trispectrum Values = '',i9/)')ot
c-------
      jf=iu-il+1
      ntr=nit+ot
      ns=n-1
      np=p-1
      nq=max(p,ntr)-1
      na=jf-1
      nn=ntr-1
      q=lb-1
c  Allocate space
      allocate(a(0:p),b(0:q),c(0:np),f(0:lb-1),in(0:nq),r(0:ns),rt(0:np)
     *,s(0:np),tr(0:nn),ti(0:nn),t1(0:nn),t2(0:nn),t3(0:nn),sp(0:na),v(0
     *:4*lb+14),w(0:ns),u(0:na,0:na,0:k3u),x(0:ns,0:nc-1),z(0:na,0:na,0:
     *k3d),stat=ibad)
      if(ibad/=0) then
       write(3,*)'Unable to allocate work space for arrays in top'
       write(6,*)'Unable to allocate work space for arrays in top'
       write(9,*)'Unable to allocate work space for arrays in top'
       stop
      endif
      call tri(a,b,c,f,in,r,rt,s,tr,ti,t1,t2,t3,sp,v,w,u,x,z,id,is,ie,n,
     *lb,il,iu,jf,iunit,kcol,k3u,k3d,nc,ntr,p,sr,pt,e,res,pv,qv,su,fmt,r
     *b,npt,nnt,ot,3)
c  Start time
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,gt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      deallocate(a,b,c,f,in,r,rt,s,sp,tr,ti,t1,t2,t3,u,v,w,x,z)
      return
c---------
  1   write(6,*)'**** Error on opening data file ',io1
      write(6,*)fmt
      write(9,*)'**** Error on opening data file ',io1
      write(9,*)fmt
      return
  2   write(6,*)'**** Error on opening ',fname,' ',io2
      write(9,*)'**** Error on opening ',fname,' ',io2
      return
  3   write(6,*)'**** Error on opening ',fname,' ',io2
      write(6,*)'**** Error on opening ',fname,' ',io2
      write(9,*)'**** Error on opening ',fname,' ',io2
      write(9,*)'**** Error on opening ',fname,' ',io2
  4   write(6,*)'**** Error on opening out file ',fname,' ',io4
      write(9,*)'**** Error on opening out file ',fname,' ',io4
      return
      end
c
c===========================================
c
      subroutine tri(a,b,c,f,in,r,rt,s,tr,ti,t1,t2,t3,sp,v,w,u,x,z,id,is
     *,ie,n,lb,il,iu,jf,iunit,kcol,k3u,k3d,nc,ntr,p,sr,pt,e,res,pv,qv,su
     *,fmt,rb,npt,nnt,ot,io)
c**********************************************************************
c  Input: x - data array, su - sampling unit, sr - sampling rate
c   fl - bin no. of lower freq, fu - bin no. of upper freq of passband
c   lb - bandwidth is 1/lb, res - resamples, pt - % taper, qv - test quantile,
c   pv - p value threshold, nit=npt (k3>0)+nnt (k3<0) - IT size, ot - OT size,
c   ntr=nit+ot - PD, nc - data columns in data file, kcol - column used
c======================================================================
      allocatable::ts
      real*8 b(lb),x(n,nc),r(n),s(p),u(jf,jf,k3u+1),w(n),z(jf,jf,k3d+1),
     *c(p),sp(jf),tr(ntr),ti(ntr),v(4*lb+15),am,sd,c3,c4,max,min,ax,sg
      real ft(1),f0(1),frv(1),fta(1),qv,tt,t0,trv,tal
      complex a(p+1),rt(p)
      complex*16 f(lb)
      integer t1(ntr),t2(ntr),t3(ntr),in(ntr),ot,p,po,q,res
      character*80 id
      character*50 fmt
      character*25 ts(:)
      character*10 su
      character mk
      intent(in) is,ie,n,lb,il,iu,jf,iunit,kcol,k3u,k3d,nc,ntr,p,sr,pt,e
     *,res,pv,qv,su,id,fmt,rb,npt,nnt,ot,io
c****************************
c  Parse format line in buf to see if there is a 'a' for time marks
      mk='a'
      ia=index(fmt,'a')
      if(ia>0) then
       mk='m'
      endif
c-------
      nmk=1
      if(mk=='m') then
       nmk=n
      endif
      allocate(ts(0:nmk-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for ts in tri'
       write(6,*)'Unable to allocate work space for ts in tri'
       stop
      endif
c  Initialize v for lb
      call cffti(lb,v)
c  Write header for summary stats in file 3
      write(io,'(67(''*''))')
      call wr2(id,su,sr,e,n,lb,pt,qv,is,ie,il,iu,res,rb,iunit,io)
      write(io,'('' Format: '',a50)') fmt
c-------
      call datc(x,ts,fmt,mk,nmk,is,ie,nc,2,io)
c  Compute statistics for kcol data
      call dstat(n,x(1,kcol),ax,sg,c3,c4,max,min,io)
      write(io,'(60(''=''))')
      write(io,'(10x,''Descriptive Statistics of Data'')')
      call dwrstat(io,ax,sg,c3,c4,max,min)
      po=0
      if(p>0) then
c  Least squares AR(p) - residuals in array r, coefficients in c
       call dfar(x(1,kcol),c,s,r,in,r2,sig,n,p,po,pv,ier,io)
       if(po>0) then
c  Compute statistics of n-po residuals
        ip=po+1
        call dstat(n-po,r(ip),am,sd,c3,c4,max,min,io)
c  Standardize residuals. First po = 0 from far
        do k=1,n
         r(k)=(r(k)-am)/sd
        enddo
c  White test
        call dport(r,n,e,ct,io)
        write(io,'(60(''=''))')
        write(io,'(10x,''Descriptive Statistics of Residuals'')')
        call dwrstat(io,am,sd,c3,c4,max,min)
        write(io,'(/7x,''White Noise Test Statistic p-value = '',f7.3/)'
     *)ct
        write(io,'(60(''*''))')
        call outar(c,s,a,rt,in,r2,sig,su,sr,pv,p,po,iunit,io)
c  Trispectrum of r
        call trisrun(r,b,w,f,in,sp,tr,ti,t1,t2,t3,u,v,z,n,lb,il,iu,jf,pt
     *,qv,tt,t0,trv,tal,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,io)
        if(res>0) then
c  Resample r where v initialized for fft with lb
         call resampl(r,b,f,in,sp,tr,ti,t1,t2,t3,u,v,w,z,n,lb,il,iu,jf,p
     *t,qv,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,ft,f0,frv,fta,tt,t
     *0,trv,tal,io)
        endif
         elseif(po==0) then
c  Standardize x => r & resample r
          do k=1,n
           r(k)=(x(k,kcol)-ax)/sg
          enddo
          call trisrun(r,b,w,f,in,sp,tr,ti,t1,t2,t3,u,v,z,n,lb,il,iu,jf,
     *pt,qv,tt,t0,trv,tal,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,io)
          if(res>0) then
           call resampl(r,b,f,in,sp,tr,ti,t1,t2,t3,u,v,w,z,n,lb,il,iu,jf
     *,pt,qv,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,ft,f0,frv,fta,tt
     *,t0,trv,tal,io)
          endif
           else
            write(io,*)'p in AR(p) call is incorrect ',po
            write(6,*)'p in AR(p) call is incorrect ',po
            write(9,*)'p in AR(p) call is incorrect ',po
            stop
       endif
c  End of branch on po in AR fit
      endif
c  Trispectrum of x
      if(p==0) then
       call trisrun(x(1,kcol),b,w,f,in,sp,tr,ti,t1,t2,t3,u,v,z,n,lb,il,i
     *u,jf,pt,qv,tt,t0,trv,tal,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3
     *d,io)
      endif
c-------
      nit=npt+nnt
      call outtris(t1,t2,t3,tr,ti,in,ft,f0,frv,fta,nit,su,rb,alm,tt,t0,t
     *rv,tal,res,iunit,io)
c  Summed on one frequency plot file
      call wrttris(t1,t2,tr,u,il,jf,nit,k3u+1,rb,iunit,7)
      call wrttris(t1,t2,tr,z,il,jf,nit,k3d+1,rb,iunit,4)
      call flush(io)
      call flush(4)
      call flush(7)
c-------
      deallocate(ts)
      return
      end
c
c===========================================
c
      subroutine trisrun(x,b,w,f,in,sp,tr,ti,t1,t2,t3,u,v,z,n,lb,il,iu,j
     *f,pt,qv,tt,t0,trv,tal,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,i
     *o)
c********************************************
c  Input: x - data, n - sample size, lb - frame size
c  pt - taper, il - low freq bin of passband, iu - high freq bin of passband
c  Output: Phases and amplitudes of trispectrum (tr,ti)
c  Counters t1,t2,t3 for trispectrum, alm - average noncentrality lamda
c  tt- CDF value for the qv nonlinearity test, t0 - CDF value for the zero
c  trispectrum test, trv - CDF value for the time reversibility test
c  tal - aliasing test CDf
c  nit=npt (k3>0)+nnt (k3<0) - IT size, ot - OT size, ntr=nit+ot - PD
c============================================================
      allocatable::y
      real*8 x(n),w(n),u(jf,jf,k3u+1),z(jf,jf,k3d+1),b(lb),y(:),sp(jf),t
     *r(ntr),ti(ntr),v(4*lb+15),fro(1),zv
      real fri(1),tt,t0,trv,tal,alm,qv
      complex*16 f(0:lb-1)
      integer t1(ntr),t2(ntr),t3(ntr),in(ntr),ot,res
      character*10 su
      intent(in) x,n,lb,il,iu,jf,pt,qv,v,rb,iunit,res,su,npt,nnt,ot,ntr,
     *k3u,k3d,io
      intent(out) alm,tt,t0,trv,tal,u,z
c********************************************
c  No. in IT
      nit=npt+nnt
      allocate(y(0:nit-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for array y'
       write(6,*)'Unable to allocate work space for array y'
       stop
      endif
c  Set indices for t1,t2,t3
      call t123(t1,t2,t3,il,iu,lb,ntr,io)
      call triloop(x,b,tr,ti,t1,t2,t3,sp,f,u,v,w,z,n,lb,il,iu,jf,npt,nnt
     *,ot,ntr,pt,t0,trv,tal,alm,su,k3u,k3d,iunit,res,io)
      do k=1,nit
       y(k-1)=1.d0-tr(k)
      enddo
c  Sort p - values in tr => u
      fri(1)=qv
      call dsort(in,y,nit,1,fri,fro,'B',ier)
      zv=(fro(1)-dble(qv))/dsqrt(dble(qv)*(1.d0-dble(qv))/dfloat(nit))
c  Nonlinearity qv quantile test cdf
      tt=dcdg(zv)
      deallocate(y)
c-------
      return
      end
c
c===========================================
c
      subroutine triloop(x,a,tr,ti,t1,t2,t3,sp,f,u,v,w,z,n,lb,il,iu,jf,n
     *pt,nnt,ot,ntr,pt,t0,trv,tal,alm,su,k3u,k3d,iunit,res,io)
c********************************************
c  Input: x - data window, n - sample size, lb - frame size, ntr - trispecs
c   pt - taper, il - low freq bin of passband, iu - high freq bin of passband
c   v - fft initialization, counters t1,t2,t3 for trispectrum, nnt - k3 < 0
c  Output: Phases and amplitudes of trispectrum (tr,ti), alm - noncentrality
c   t0 - CDF value for the trispec=0 test, trv - CDF value for the time
c   reversibility test, tal - aliasing test
c============================================================
      real*8 a(lb),x(n),w(n),u(jf,jf,k3u+1),z(jf,jf,k3d+1),sp(jf),tr(ntr
     *),ti(ntr),v(4*lb+15),sw
      real t0,trv,tal,alm
      complex*16 f(0:lb-1)
      character*10 su
      integer t1(ntr),t2(ntr),t3(ntr),ot,res
      intent(in) x,v,t1,t2,t3,n,lb,il,iu,jf,npt,nnt,ot,ntr,pt,su,k3u,k3d
     *,iunit,res,io
      intent(out) alm,tr,ti,sp,t0,trv,tal,u,z
c********************************************
c  Initialize taper
      nt=0
      if(pt>0.0) then
       call dtaper(w,lb,pt,sw,nt,io)
c  If no. of tapered values at ends < 3 set sw=1
       if(nt<3) then
        sw=1.d0
       endif
        else
         nt=0
         sw=1.d0
      endif
c-------
      call p234(x,sp,f,v,w,a,tr,ti,n,lb,ntr,ot,il,iu,jf,nt,io)
c  Normalized trispectrum probs
      call triscale(tr,ti,t1,t2,t3,u,z,sp,n,lb,il,jf,ntr,npt,nnt,t0,trv,
     *tal,alm,sw,k3u,k3d,nt,io)
c  Add over f3
      call scalu(u,z,il,iu,jf,k3u,k3d,io)
c------
      return
      end
c
c===========================================
c
      subroutine p234(x,sp,f,v,w,a,tr,ti,n,lb,ntr,ot,il,iu,jf,nt,io)
c********************************************************
c  Periodogram & trispectrum frame loop
c   t1, t2, t3 - IT bins, nit=npt (k3>0)+nnt (k3<0) - IT size, ot - OT size,
c   nfrm - no. of frames, (io,iu) band counters, nt - taper
c========================================================
      real*8 x(n),a(lb),w(lb),sp(jf),tr(ntr),ti(ntr),v(4*lb+15),dlb,dnu
      complex*16 f(0:lb-1),c4,cnu
      integer t1(ntr),t2(ntr),t3(ntr),sr,ot
      intent(in) x,v,w,lb,n,ntr,ot,il,iu,jf,io
      intent(out) sp,tr,ti
c********************************************
c  Initialize arrays
      do k=1,jf
       sp(k)=0.d0
       tr(k)=0.d0
       ti(k)=0.d0
      enddo
      do k=jf+1,ntr
       tr(k)=0.d0
       ti(k)=0.d0
      enddo
c-------
      dnu=dfloat(n)
c  Complex n
      cnu=dcmplx(dnu)
      nfrm=n/lb
      ks=il-1
c  Start loop over frames
      do kf=1,nfrm
       kp=(kf-1)*lb
       if(nt>=3) then
c  Taper frame kb
        do k=1,lb
c  x pointer
         ka=kp+k
         a(k)=w(k)*x(ka)
        enddo
         else
c  Shift to a
          do k=1,lb
           ka=kp+k
           a(k)=x(ka)
          enddo
       endif
c  w frame => f
       do j=1,lb
        f(j-1)=dcmplx(dble(a(j)))
       enddo
c  FFT of frame
       call cfftf(lb,f,v)
c-------
       sr=0
       do k1=il,iu
c  Shift to counter j1 in sp
        j1=k1-ks
c  |f|**2/nu
        sp(j1)=sp(j1)+dreal(f(k1)/cnu*f(lb-k1))
c  k3 > 0
        do k2=il,k1
         if((k1+k2)<iu) then
          do k3=il,min(k2,iu-k1-k2)
           k4=-k1-k2-k3
           sr=sr+1
           c4=f(k1)/cnu*f(k2)*f(k3)*f(lb+k4)
           tr(sr)=tr(sr)+dreal(c4)
           ti(sr)=ti(sr)+dimag(c4)
          enddo
         endif
        enddo
       enddo
c  k3 < 0
       do k1=il,iu
        do k2=il,k1
         mk=max(il,k1+k2-iu)
         do k3=-mk,-k2,-1
          k4=-k1-k2-k3
c  Skip submanifolds
          if(((k3+k4)/=0).and.((k1+k3)/=0.and.(k2+k4)/=0).and.((k1+k4)/=
     *0.and.(k2+k3)/=0)) then
           sr=sr+1
           c4=f(k1)/cnu*f(k2)*f(lb+k3)*f(lb+k4)
           tr(sr)=tr(sr)+dreal(c4)
           ti(sr)=ti(sr)+dimag(c4)
          endif
         enddo
        enddo
       enddo
c-------
       kot=sr
c  OT k3 > 0
       do k1=il,min(lb/6,iu)
        do k2=il,k1
         do k3=il,min(k2,iu-k1-k2)
          nkot=k1+k2+k3
          kot=kot+1
          c4=f(k1)/cnu*f(k2)*f(k3)*f(nkot)
          tr(kot)=tr(kot)+dreal(c4)
          ti(kot)=ti(kot)+dimag(c4)
         enddo
        enddo
       enddo
c  OT for k3 < 0
       do k1=max(il,lb/4),iu
        do k2=il,k1
         mk=max(il,k1+k2-iu)
         do k3=-mk,-k2,-1
          if((k2+k3)/=0) then
           nkot=k1+k2+k3
           kot=kot+1
           c4=f(k1)/cnu*f(k2)*f(lb+k3)*f(nkot)
           tr(kot)=tr(kot)+dreal(c4)
           ti(kot)=ti(kot)+dimag(c4)
          endif
         enddo
        enddo
       enddo
       if(sr/=ntr-ot) then
        write(io,*) 'IT count ',ntr-ot,' /= sr ',sr,' in p234'
        stop
       endif
       if(kot/=ntr) then
        write(io,*) 'Array count ',kot,' /= ntr ',ntr,' in p234'
        stop
       endif
      enddo
c  End of loop on frames
c-------
      return
      end
c
c===========================================
c
      subroutine triscale(tr,ti,t1,t2,t3,u,z,sp,n,lb,il,jf,ntr,npt,
     *nnt,t0,trv,tal,alm,sw,k3u,k3d,nt,io)
c****************************************************************************
c     Computes trispectrum  12-21-2007
c     t1, t2, t3 - IT bins, ntr & not - sizes of tr & ti arrays for IT & OT
c     t0 - CDF value for the zero trispectrum null, nfrm - frames, nnt - k3 < 0
c     trv - CDF value for the zero imaginary part
c     alm - average noncentrality lamda of the trispectral estimate
c============================================================================
      real rdb,t0,trv,tal,alm
      real*8 u(jf,jf,k3u+1),z(jf,jf,k3d+1),sp(jf),tr(ntr),ti(ntr),ar,ai,
     *a1,a2,as,beta,h,scale,sw
      integer t1(ntr),t2(ntr),t3(ntr),n,lb,il,jf,ntr,npt,nnt,nt,ks,nt0,n
     *as,nit,k,k1,k2,k3,k4,ndb,ot,k3d,k3u,iuu,io
      intent(in) t1,t2,t3,sp,n,lb,il,jf,ntr,npt,nnt,sw,k3u,k3d,nt,io
      intent(out) u,z,t0,trv,tal,alm
      intent(in out) tr,ti
c********************************************
      scale=dsqrt(dfloat(n/lb)*2.d0)/dfloat(lb)
      ks=il-1
      iuu=jf+il-1
c  Initialize count for spec=0 values
      nt0=0
      nas=0
      a1=0.d0
      a2=0.d0
      nit=npt+nnt
c  k3>0
      do k=1,npt
       k1=t1(k)
       k2=t2(k)
       k3=t3(k)
       k4=-(k1+k2+k3)
       call scaledet(k1,k2,k3,k4,iuu,beta,io)
       h=beta*sp(k1-ks)*sp(k2-ks)*sp(k3-ks)*sp(-k4-ks)
       if(nt>=3) then
c  Taper scale
        h=h/(sw*sw*sw*sw)
       endif
       if(h>10**(-25)) then
        tr(k)=tr(k)/dsqrt(h)*scale
        ti(k)=ti(k)/dsqrt(h)*scale
c  ReTris**2
        ar=tr(k)*tr(k)
c  ImTris**2
        ai=ti(k)*ti(k)
c  Sum of squares
        as=ar+ai
        a1=a1+ar
        a2=a2+ai
c  Phases of Tris
        ti(k)=datan2(ti(k),tr(k))*(57.2958)
c  |Tris|**2
        tr(k)=as
c  u trispectrum
        u(t1(k)-ks,t2(k)-ks,t3(k)-ks)=as
         else
          nt0=nt0+1
          tr(k)=0.d0
          ti(k)=0.d0
          u(t1(k)-ks,t2(k)-ks,t3(k)-ks)=0.0
       endif
      enddo
c  k3<0
      do k=npt+1,nit
       k1=t1(k)
       k2=t2(k)
       k3=t3(k)
       k4=-(k1+k2+k3)
       call scaledet(k1,k2,k3,k4,iuu,beta,io)
       h=beta*sp(k1-ks)*sp(k2-ks)*sp(-k3-ks)*sp(-k4-ks)
       if(nt>=3) then
c  Taper scale
        h=h/(sw*sw*sw*sw)
       endif
       if(h>10**(-25)) then
        tr(k)=tr(k)/dsqrt(h)*scale
        ti(k)=ti(k)/dsqrt(h)*scale
c  ReTris**2
        ar=tr(k)*tr(k)
c  ImTris**2
        ai=ti(k)*ti(k)
c  Sum of squares
        as=ar+ai
        a1=a1+ar
        a2=a2+ai
c  Phases of Tris
        ti(k)=datan2(ti(k),tr(k))*(57.2958)
c  |Tris|**2
        tr(k)=as
c  z trispectrum
        z(t1(k)-ks,t2(k)-ks,-t3(k)-ks)=as
         else
          nt0=nt0+1
          tr(k)=0.d0
          ti(k)=0.d0
          z(t1(k)-ks,t2(k)-ks,-t3(k)-ks)=0.0
        endif
      enddo
c-------
c  No of terms with non-zero spectral products
      ndb=nit-nt0
      rdb=float(ndb)
      if(ndb<12) then
       write(io,'(40(''*'')/)')
       write(io,'('' Valid trispectrum values = '',i1,'' < 12'')')ndb
       do k=1,nit
        tr(k)=0.d0
       enddo
       return
      endif
c  Noncentrality lamda
      alm=max(sngl(a1+a2)/rdb-2.,0.0)
c  Uniform cdf values under linear null hypothesis
      do k=1,nit
       tr(k)=p2chi(sngl(tr(k)),alm,io)
      enddo
c  CDF statistic U(0,1) of normalized a1+a2 - zero trispectrum test
      call cdchi(sngl(a1+a2),2.*rdb,t0,io)
c  CDF statistic U(0,1) of normalized a2 - time reversibility test
      call cdchi(sngl(a2),rdb,trv,io)
c  Aliasing test
      ot=ntr-nit
      a1=0.d0
      a2=0.d0
      do k=nit+1,ntr
       k1=t1(k)
       k2=t2(k)
       k3=t3(k)
       k4=lb-(k1+k2+k3)
       call scaledet(k1,k2,k3,k4,iuu,beta,io)
       h=beta*sp(k1-ks)*sp(k2-ks)*sp(abs(k3)-ks)*sp((lb-k4)-ks)
       if(h>10**(-25)) then
        tr(k)=tr(k)/dsqrt(h)*scale
        ti(k)=ti(k)/dsqrt(h)*scale
c  ReTris**2
        ar=tr(k)*tr(k)
c  ImTris**2
        ai=ti(k)*ti(k)
c  Sum of squares
        a1=a1+ar
        a2=a2+ai
c  |Tris|**2
        tr(k)=ar+ai
         else
          nas=nas+1
       endif
      enddo
c  CDF statistic U(0,1) of aliasing test
      call cdchi(sngl(a1+a2),2.*float(ot-nas),tal,io)
c-------
      return
      end
c
c===========================================
c
      subroutine scaledet(k1,k2,k3,k4,iu,beta,io)
c************************************************
c     Determine scale factos given frequency index values
c     store value as beta
c     See Table II - Dalle Molle and Hinich (1995) - JASOA p. 2970
c================================================
      real*8 beta
      integer k1,k2,k3,k4,iu,io
      intent(in) k1,k2,k3,k4,iu,io
      intent(out) beta
c***********************************
c  Default case
      beta=1.d0
c  case 1
      if((k1==k2).and.(k2/=k3).and.(k3/=k4)) then
       beta=2.d0
      endif
c     case 2
      if((k1==k2).and.(k2/=k3).and.(k3==k4)) then
       beta=4.d0
      endif
c  case 3
      if((k1==k2).and.(k2==k3).and.(k3/=k4)) then
       beta=6.d0
      endif
c  case 4
      if((k1==k2).and.(k2==k3).and.(k3==k4)) then
       beta=24.d0
      endif
c  case 5
      if((k1==k2).and.(k2==(-k3)).and.((-k3)==(-k4))) then
       beta=20.d0
      endif
c  case 6
      if((k1==k2).and.(k2==k3).and.(k3/=k4).and.(k4==iu)) then
       beta=6.d0
      endif
c  case 7
      if((k1==k2).and.(k2==k3).and.(k3==k4).and.(k4==iu)) then
       beta=96.d0
      endif
      return
      end
c
c===========================================
c
      subroutine scalu(u,z,il,iu,jf,k3u,k3d,io)
c*******************************************
c  Set u & z for writtris
c===========================================
      real*8 u(jf,jf,k3u+1),z(jf,jf,k3d+1),stor
      intent(in) il,iu,jf,k3u,k3d
      intent(in out) u,z
c*******************************************
      n3u=k3u+1
      n3d=k3d+1
      ks=il-1
c  Zero k2>k1 part of u & z
      do j1=1,jf
       do j2=j1+1,jf
        u(j1,j2,n3u)=0.d0
        z(j1,j2,n3d)=0.d0
       enddo
      enddo
c  k3 > 0
      do k1=il,iu
       do k2=il,k1
c  Sum over k3
        stor=0.d0
        ns=0
        if((k1+k2)<iu) then
         do k3=il,min(k2,iu-k1-k2)
          stor=stor+dble(u(k1-ks,k2-ks,k3-ks))
          ns=ns+1
         enddo
c  Chi(2*k3u) => u
         if(ns>0) then
          call cdchi(sngl(stor),2.0*float(ns),uh,io)
           else
            uh=0.0
         endif
         u(k1-ks,k2-ks,n3u)=dble(uh)
          else
           u(k1-ks,k2-ks,n3u)=0.d0
        endif
       enddo
      enddo
c  k3 < 0
      do k1=il,iu
       do k2=il,k1
        mk=max(il,k1+k2-iu)
        stor=0.d0
        ns=0
        do k3=-mk,-k2,-1
         k4=-k1-k2-k3
c  Skip submanifolds
         if(((k3+k4)/=0).and.((k1+k3)/=0.and.(k2+k4)/=0).and.((k1+k4)/=0
     *.and.(k2+k3)/=0)) then
          stor=stor+dble(z(k1-ks,k2-ks,-k3-ks))
          ns=ns+1
         endif
        enddo
c  Chi(2*k3d) =>z(:,:,1)
        if(ns>0) then
         call cdchi(sngl(stor),2.*float(ns),uh,io)
          else
           uh=0.0
        endif
        z(k1-ks,k2-ks,n3d)=dble(uh)
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine t123(t1,t2,t3,fl,fu,lb,ntr,io)
c*******************************************
c  Set indices for t1, t2, t3
c===========================================
      integer t1(ntr),t2(ntr),t3(ntr),fl,fu,sr
      intent(in) fl,fu,lb,ntr
      intent(out) t1,t2,t3
c*******************************************
      do k=1,ntr
       t1(k)=0
       t2(k)=0
      enddo
      sr=0
c  k3 > 0
      do k1=fl,fu
       do k2=fl,k1
        if((k1+k2)<fu) then
         do k3=fl,min(k2,fu-k1-k2)
          sr=sr+1
c  Set indices
          t1(sr)=k1
          t2(sr)=k2
          t3(sr)=k3
         enddo
        endif
       enddo
      enddo
c  k3 < 0
      do k1=fl,fu
       do k2=fl,k1
        mk=max(fl,k1+k2-fu)
        do k3=-mk,-k2,-1
         k4=-k1-k2-k3
c  Skip submanifolds
         if(((k3+k4)/=0).and.((k1+k3)/=0.and.(k2+k4)/=0).and.((k1+k4)/=0
     *.and.(k2+k3)/=0)) then
c  Set indices
          sr=sr+1
          t1(sr)=k1
          t2(sr)=k2
          t3(sr)=k3
         endif
        enddo
       enddo
      enddo
c-------
c  Start counting - OT
      ot=sr
c  OT k3 > 0
      do k1=fl,min(lb/6,fu)
       do k2=fl,k1
        do k3=fl,min(k2,fu-k1-k2)
         ot=ot+1
         t1(ot)=k1
         t2(ot)=k2
         t3(ot)=k3
        enddo
       enddo
      enddo
c  OT for k3 < 0
      do k1=max(fl,lb/4),fu
       do k2=fl,k1
        mk=max(fl,k1+k2-fu)
        do k3=-mk,-k2,-1
         if((k2+k3)/=0) then
          ot=ot+1
          t1(ot)=k1
          t2(ot)=k2
          t3(ot)=k3
         endif
        enddo
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine size(fl,fu,lb,npt,nnt,ot,k3u,k3d,io)
c*******************************************
c  Computes trispectrum loops for array size
c   npt - IT for k3>0, nnt - IT for k3<0, ot- for OT, k3u - k3+, k3d - k3-
c===========================================
      integer fl,fu,ot
      intent(in) fl,fu,lb
      intent(out) npt,nnt,ot,k3u,k3d
c*******************************************
      npt=0
      k3u=0
      k3d=0
c  k3 > 0
      do k1=fl,fu
       do k2=fl,k1
        if((k1+k2)<fu) then
         do k3=fl,min(k2,fu-k1-k2)
          npt=npt+1
          if(k3>k3u) then
           k3u=k3
          endif
         enddo
        endif
       enddo
      enddo
      nnt=0
c  k3 < 0
      do k1=fl,fu
       do k2=fl,k1
        mk=max(fl,k1+k2-fu)
        do k3=-mk,-k2,-1
         k4=-k1-k2-k3
c  Skip submanifolds
         if(((k3+k4)/=0).and.((k1+k3)/=0.and.(k2+k4)/=0).and.((k1+k4)/=0
     *.and.(k2+k3)/=0)) then
          if(-k3>k3d) then
           k3d=-k3
          endif
          nnt=nnt+1
         endif
        enddo
       enddo
      enddo
c-------
c  Start counting - OT
      ot=0
c  OT k3 > 0
      do k1=fl,min(lb/6,fu)
       do k2=fl,k1
        do k3=fl,min(k2,fu-k1-k2)
         ot=ot+1
        enddo
       enddo
      enddo
c  OT for k3 < 0
      do k1=max(fl,lb/4),fu
       do k2=fl,k1
        mk=max(fl,k1+k2-fu)
        do k3=-mk,-k2,-1
         if((k2+k3)/=0) then
          ot=ot+1
         endif
        enddo
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine resampl(x,b,f,in,sp,tr,ti,t1,t2,t3,u,v,w,z,n,lb,il,iu,j
     *f,pt,qv,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,ft,f0,frv,fta,t
     *t,t0,trv,tal,io)
c********************************************
c  Computes fi p-value thresholds for resamples=res using shuffle bootstrap
c  x - data, w - shuffle
c  ft - nonlinearity test, f0 - zero trispectrum, frv - time reversibility
c  fta - aliasing thresholds
c  tt- CDF value for the qv nonlinearity test, t0 - CDF value for the zero
c  trispectrum test, trv - CDF value for the time reversibility test
c  tal - aliasing test CDf
c============================================
      allocatable::at,a0,arv,ata,hold
      real*8 x(n),w(n),u(jf,jf,k3u+1),z(jf,jf,k3d+1),b(lb),sp(jf),tr(ntr
     *),ti(ntr),v(4*lb+15),at(:),a0(:),arv(:),ata(:)
      real ft(1),f0(1),frv(1),fi(1),fta(1),rt,r0,rrv,ral,alm,qv,tt,t0,tr
     *v,tal
      complex*16 f(0:lb-1)
      integer t1(ntr),t2(ntr),t3(ntr),in(ntr),hold(:),ot,res
      character*10 su
      intent(in) x,n,lb,il,iu,jf,pt,qv,v,rb,iunit,res,su,npt,nnt,ot,ntr,
     *k3u,k3d,tt,t0,trv,tal,io
      intent(out) ft,f0,frv,fta
c********************************************
c  Allocate space
      ns=res-1
      allocate(at(0:ns),a0(0:ns),arv(0:ns),ata(0:ns),hold(0:n-1),stat=ib
     *ad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space in resampl'
       stop
      endif
c-------
      n2=lb/2
      do kr=1,res
       k=kr-1
       call dshuf(x,w,hold,n,io)
c  Trispectrum of w
       call trisrun(w,b,x,f,in,sp,tr,ti,t1,t2,t3,u,v,z,n,lb,il,iu,jf,pt,
     *qv,rt,r0,rrv,ral,alm,iunit,res,rb,su,npt,nnt,ot,ntr,k3u,k3d,io)
       at(k)=1.d0-dble(rt)
       a0(k)=1.d0-dble(r0)
       arv(k)=1.d0-dble(rrv)
       ata(k)=1.d0-dble(ral)
c  End of resample loop
      enddo
c  Thresholds for p-values for nonlinearity test
      fi(1)=1.0-tt
      call dsort(hold,at,res,7,fi,ft,'B',ier)
c  Thresholds for p-values for zero trispectrum test
      fi(1)=1.0-t0
      call dsort(hold,a0,res,7,fi,f0,'B',ier)
c  Thresholds for p-values for time reversibility test
      fi(1)=1.0-trv
      call dsort(hold,arv,res,7,fi,frv,'B',ier)
c  Thresholds for p-values for aliasing test
      fi(1)=1.0-tal
      call dsort(hold,ata,res,7,fi,fta,'B',ier)
c-------
      deallocate(at,a0,arv,ata,hold)
      return
      end
c
c===========================================
c
      subroutine datc(x,ts,fmt,mk,nmk,is,ie,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*8 x(ie-is+1,nc)
      character*50 fmt
      character*25 ts(nmk)
      character mk
      intent(in) fmt,mk,nmk,is,ie,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(9,'(7x,''/Data & time format is incorrect''/)')
            write(9,'(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-is+1),(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-is
     *+1)
         endif
        enddo
       endif
      endif
      if(mk/='m') then
       if(is==1) then
        do i=1,ie
         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,*)'t = ',i-is+1
      write(6,*)'**** eof in read from data file '
      write(6,*)'t = ',i-is+1
      write(9,*)'**** eof in read from data file '
      write(9,*)'t = ',i-is+1
      stop
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(9,*)'** Error in read of data file ',ko
      write(9,*)'Check format',fmt
      stop
      end
c
c===========================================
c
      subroutine wr2(id,su,sr,e,n,lb,pt,qv,is,ie,il,iu,res,rb,iunit,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      real pt,qv
      integer res
      character*79 id
      character*10 su
      intent(in) id,su,sr,e,n,lb,pt,qv,is,ie,il,iu,res,rb,iunit,io
c*****************************************
      jf=iu-il+1
      rl=float(lb)
      write(io,'(1x,a79)')id
      write(io,'(/7x,''First datum read'',i7,2x,''Last datum read'',i8)'
     *)is,ie
      if(iunit>0) then
       write(io,'(/'' Sample size = '',i8,2x,''No. of frequencies in the
     * band = '',i6,2x,''No. frames = '',i6/)')n,jf,n/lb
      endif
      if(iunit==0) then
       write(io,'(/'' Sample size = '',i8,2x,''No. of periods in the ban
     *d = '',i6,2x,''No. frames = '',i6/)')n,jf,n/lb
      endif
c-------
      if(iunit==1) then
c  Lower and upper frequencies for band
       a=float(il)*rb*1.e-3
       b=float(iu)*rb*1.e-3
       write(io,'('' Sampling rate ='',f12.4,'' kHz '',4x,''Frame size '
     *',i5/)') sr,lb
       write(io,'('' Resolution Bandwidth :'',g14.4,'' Hz'',g14.4,'' mse
     *c''/)')rb,1.0/rb
       write(io,'('' Passband ('',f9.4,2x,f10.4,'' ) kHz''/)')a,b
      endif
      if(iunit==2) then
c  Lower and upper frequencies for band
       a=float(il)*rb
       b=float(iu)*rb
       write(io,'('' Sampling rate ='',f12.4,'' Hz '',4x,''Frame size ''
     *,i5/)')sr,lb
       write(io,'('' Resolution Bandwidth :'',g14.4,'' Hz'',g14.4,'' mse
     *c''/)') rb,1.0/rb
       write(io,'('' Passband ('',f9.4,2x,f10.4,'' ) Hz''/)')a,b
      endif
      if(iunit==0) then
c  Time unit
         a=rb/float(il)
         b=rb/float(iu)
         write(io,'('' Sampling interval = '',f12.4,1x,a10,2x,''Frame si
     *ze '',i5/)')sr,su,lb
         write(io,'('' Resolution Bandwidth :'',g14.4,1x,a10/)')rb,su
         write(io,'('' Passband ('',f10.2,2x,f7.4,'' )'',1x,a10/)')a,b,
     *su
      endif
c-------
c  Number of values tapered from each end
      nt=nint(pt*rl/200.)
c  If no. of tapered values at ends < 2 set nt=0
      if(nt<2) then
       nt=0
      endif
      if(pt>0.0 .and. nt>0) then
       write(io,'('' % taper = '',f5.2,4x,''No. tapered at each end = ''
     *,i5/)')pt,nt
      endif
      if(nt==0) then
       write(io,'('' % Taper yields no. of tapered values < 2''/)')
      endif
      if(res>0) then
       write(io,'(4x,''Residuals are Bootstrapped using '',i5,'' resampl
     *es''/)')res
      endif
      write(io,'(7x,''Nonlinearity Test Computed for the '',f5.2,'' Quan
     *tile''/)')qv
      write(io,'(7x,''Bandwidth Frame Exponent = '',f7.4/)')e
c-------
      return
      end
c
c===================================================
c
      subroutine dwrstat(io,am,sig,sk,c4,max,min)
c****************************************************
c  Data statistics
c============================================
      implicit none
      real*8 am,sig,sk,c4,max,min
      integer io
      intent(in) am,sig,sk,c4,max,min,io
c********************************************
      write(io,'(60(''='')/)')
      write(io,'('' Mean = '',f15.3,7x,''Std Dev = '',f15.3/)')am,sig
      write(io,'('' Skew = '',g10.4,4x,''Kurtosis = '',g10.4/)')sk,c4
      write(io,'('' Max value = '',f15.3,7x,''Min value = '',f15.3)')max
     *,min
      write(io,'(60(''+''))')
c-------
      return
      end
c
c=============================================
c
      subroutine wrttris(t1,t2,tr,u,il,jf,nit,ku,rb,iunit,io)
c*******************************************
c  Output trispectrum for plotting 9-11-2004
c===========================================
      real*8 tr(nit),u(jf,jf,ku)
      integer t1(nit),t2(nit)
      intent(in) t1,t2,tr,u,il,jf,nit,ku,rb,iunit
c********************************************
      ks=il-1
      iu=il-1+jf
      do k2=il,iu
       if(iunit>0) then
        write(io,'(:,10000(f9.4))')(u(k1-ks,k2-ks,ku),k1=il,iu)
         else
          write(io,'(:,10000(f14.3))')(u(k1-ks,k2-ks,ku),k1=il,iu)
       endif
      enddo
      if(iunit==1) then
       write(io,'(:,10000(f9.3))')(float(k)*rb*1.e-3,k=il,iu)
        elseif(iunit==2) then
         write(io,'(:,10000(f9.3))')(float(k)*rb,k=il,iu)
          else
           write(io,'(:,10000(f14.3))')(rb/float(k),k=il,iu)
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine outtris(t1,t2,t3,tr,ti,in,ft,f0,frv,fta,nit,su,rb,alm,t
     *t,t0,trv,tal,res,iunit,io)
c*******************************************
c  Output large trispectral values  nit=npt (k3>0)+nnt (k3<0) - IT size
c===========================================
      allocatable::ak,r
      real*8 tr(nit),ti(nit),r(:),fo(11),ab,rf
      real fi(11),ft(1),f0(1),frv(1),fta(1),tt,t0,trv,tal,alm
      integer t1(nit),t2(nit),t3(nit),in(nit),ak(:),ot,res
      character*10 su
      data fi/.001,.01,.05,.1,.25,.5,.75,.9,.95,.99,.999/
      intent(in) t1,t2,t3,tr,ti,in,ft,f0,frv,fta,nit,su,rb,alm,tt,t0,trv
     *,tal,res,iunit,io
c********************************************
      allocate(ak(0:nit-1),r(0:nit-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work ak in outtris'
       stop
      endif
c  Write p-values for tests
      write(io,'(/2x,''p-value = '',f7.5,'' for the quantile trispectrum
     *based nonlinearity test'')')1.0-tt
      write(io,'(/2x,''p-value = '',f7.5,'' for the time reversibility t
     *est'')')1.0-trv
      write(io,'(/2x,''p-value = '',f7.5,'' for the zero trispectrum tes
     *t'')')1.0-t0
      write(io,'(/2x,''p-value = '',f7.5,'' for the aliasing test''/)')
     *1.0-tal
      if(res>0) then
       write(io,'(2x,70(''*'')/)')
c  Write shuffled p-values for tests
       write(io,'(2x,''Shuffled p-value = '',f7.4,'' for the quantile tr
     *ispectrum based nonlinearity test'')')ft(1)
       write(io,'(/2x,''Shuffled p-value = '',f7.4,'' for the time rever
     *sibility test'')')frv(1)
       write(io,'(/2x,''Shuffled p-value = '',f7.4,'' for the zero trisp
     *ectrum test'')')f0(1)
       write(io,'(/2x,''Shuffled p-value = '',f7.4,'' for the aliasing t
     *est''/)')fta(1)
      endif
c-------
      write(io,'(2x,70(''='')/)')
c  Write lamda for chi2 of sum of squares
      write(io,'(12x,''Mean noncentrality parameter = '',f10.2/)')alm
      do k=1,nit
       r(k-1)=sngl(tr(k))
      enddo
      call dsort(ak,r,nit,11,fi,fo,'B',ier)
c  Write fractiles
      write(io,'(11x,''Fractiles of Normalized Trispectral Probabilities
     *''/)')
      write(io,'(4x,''0.001'',4x,''0.01'',5x,''0.05'',5x,''0.10'',5x,''0
     *.25'',5x,''0.50'',5x,''0.75'',5x,''0.90'',5x,''0.95'',5x,''0.99'',
     *5x,''0.999''/)')
      write(io,'(2x,11(f7.3,2x)/)') (fo(k),k=1,11)
      do k=1,11
       rf=dble(fi(k))
       ab=dsqrt(rf*(1.d0-rf)/dfloat(nit))
       ab=abs(fo(k)-rf)/ab
       fo(k)=1.d0-dcdg(ab)
      enddo
      write(io,'(2x,11(f7.3,2x),'' p-values'')')(fo(k),k=1,11)
      write(io,'(2x,70(''='')/)')
c-------
c  Output top trispectral values
      nout=0
      pt=1.0/float(nit)
      do j=1,nit
       ji=in(j)
       if(sngl(tr(ji))<pt) then
        nout=nout+1
        ak(nout-1)=ji
       endif
      enddo
      if(nout>0) then
       write(io,'(4x,i5,'' Trispectral p-values < '',f8.5)')nout,pt
       write(io,'(2x,''% of Trispectral Probabilites = '',f6.3,2x,''Stan
     *dard Deviation = '',f6.3)')(float(nout)/float(nit))*100.0,sqrt(pt*
     *(1.0-pt)/float(nit))*100.0
       else
        return
      endif
c-------
      write(io,'(1x,60(''='')/)')
      if(iunit==1) then
       write(io,'(12x,'' Frequencies in kHz''/)')
       write(io,'(5x,''f(1)'',7x,''f(2)'',7x,''f(3)'',7x,''f(4)'',11x,''
     *1-Pr'',8x,''Phase''/)')
        elseif(iunit==2) then
         write(io,'(12x,'' Frequencies in Hz''/)')
         write(io,'(5x,''f(1)'',7x,''f(2)'',7x,''f(3)'',7x,''f(4)'',11x,
     *''1-Pr'',8x,''Phase''/)')
      endif
      if(iunit==1) then
       do j=1,nout
        ji=ak(j-1)
        s1=float(t1(ji))*rb*1.e-3
        s2=float(t2(ji))*rb*1.e-3
        s3=float(t3(ji))*rb*1.e-3
        s4=float(-t1(ji)-t2(ji)-t3(ji))*rb*1.e-3
        pr=tr(ji)
        if(s1>1.e-2.or.s2>1.e-2) then
         write(io,'(1x,4(f9.3,2x),7x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti(j
     *i)
          else
           write(io,'(1x,4(g9.3,2x),7x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti
     *(ji)
        endif
       enddo
      endif
      if(iunit==2) then
       do j=1,nout
        ji=ak(j-1)
        s1=float(t1(ji))*rb
        s2=float(t2(ji))*rb
        s3=float(t3(ji))*rb
        s4=float(-t1(ji)-t2(ji)-t3(ji))*rb
        pr=tr(ji)
        if(s1>1.e-2.or.s2>1.e-2) then
         write(io,'(1x,4(f9.3,2x),7x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti(j
     *i)
          else
           write(io,'(1x,3(g9.3,2x),7x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti
     *(ji)
        endif
       enddo
      endif
      if(iunit==0) then
c  Time unit
       write(io,'(17x,''Periods in '',a10/)')su
       write(io,'(10x,''Period 1'',8x,''Period 2'',8x,''Period 3'',8x,''
     *Period 4'',9x,''1-Pr'',8x,''Phase''/)')
       do j=1,nout
        ji=ak(j-1)
        s1=rb/float(t1(ji))
        s2=rb/float(t2(ji))
        s3=rb/float(t3(ji))
        s4=rb/float(-t1(ji)-t2(ji)-t3(ji))
        pr=tr(ji)
        if(s1>1.e-2.or.s2>1.e-2) then
         write(io,'(2x,4(f14.3,2x),7x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti(
     *ji)
         else
          write(io,'(7x,4(g11.3,5x),5x,f8.6,5x,f6.1/)')s1,s2,s3,s4,pr,ti
     *(ji)
        endif
       enddo
      endif
c--------
      deallocate(ak,r)
      return
      end