program data
c********************************************************
c  Version 3-1-2010    Melvin Hinich
c  File 1 - Rundata.cnl, File 2 - data file, File 3 - fname.out for each fname
c  File 4 - Rearranged data if called for
c
c  This program can read columns of data from a time series data file or one
c  column from a number of time series data files.
c  It can perform several transformations on the data from each column.
c  For example the raw data or the transformed data can then be filtered using
c  a bandpass filter. The filtered data can be decimated (skip sampled) at
c  twice the highest frequency of the band.
c
c  The data filtered or not can be passed to program that fits an AR(p) model
c  using an iterative fitting method. Using a probability value threshold
c  set by the user, the p-value of the t values for each lag parameter is
c  compared with the threshold and the lag is retained if and only if the
c  p-value for the lag parameter estimate is smaller than the p-value
c  threshold. Otherwise the lag is dropped and a new least squares fit is made.
c
c  If the AR fit is selected the residuals are written to the file
c  fname-co_out.data  where fname is the data file name and co is the column
c  number.
c
c  The transformations are: log, differences in logs (rates), first differences
c  detrend by a curvilinear least squares
c
c  Then the raw or transformed data can have its outliers trimmed as explained
c  below or clipped to 1/-1 around the sample mean.
c
c=================================================
c  Parameters in control file Rundata.cnl. Use spaces between entries!
c  A delimiter symbol * must be placed on each line.
c  Place comments on the lines of rundata.cnl file after the symbol *.
c  Rundata.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.data where
c  'filename' is the name of the data file used.
c===================================================================
c
c  Lines for the control file Rundata.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 control files called by rundata.cnl. The integer to the right
c  of the = sign in rundata.cnl is read by the program and sets the number
c  of runs.
c  In any *.cnl called by Rundata.cnl, !ru is the wildcard for the upper
c  bandwidth. The real number to the right of !r= is the upper bandwidth
c  selected in the control file.
c
c  ALWAYS use 1.0 rather than 1. since the program searches for the number to
c  the left of the period in the parameter.
c
c  The delimiter symbol * must be placed on each line after the parameters.
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  delimiter symbol *. These are not read by the program.
c
c===================================================================
c
c  Path & name of control file=\prog\data-file.cnl     * cnl file
c
c  The path & name of the control (cnl) file that passes parameters to the
c  program.
c
c============================
c
c  Number of data files=!nfil   (!nfil - integer wildcard)
c
c  If !nfil > 1 there has to be !nfil file names in rundata.cnl AND each data
c  file MUST have the same number of samples.
c  The maximum number of files <=40
c
c============================
c
c  Suppose that !nfil=1 & there are more than one column in the data file.
c  If you want to put the transformed data files that are selected in the
c  control file described below type Yes or yes after the = sign. For example
c
c  Output transformed columns in one file output=yes
c   OR
c  Output transformed columns in one file output=Yes
c
c  If this line is omitted in Rundata.cnl or if no or No is typed then
c  the transformed data for each column will be in a separate file.
c
c============================
c
c  The data file lines are followed by the line
c  # * End of Rundata.cnl data file reads
c
c  The program checks for the delimiter symbol # to end the reading of the
c  parameters for each control file set in Rundata.cnl
c
c============================
c
c  The loop on the k = 1,...,!nfil begins for !nfil data files.
c  Each data file name is now read.
c
c  For example if !nfile=2:
c
c  data name=\prog\data-fil1.data            * 1st data file name
c  data name=\prog\data-fil2.data            * 2nd data file name
c
c============================
c
c  The program now reads the following lines from the previously named cnl file.
c  Type the following commands for each columns followed by the line:
c  # */ End of run parameter reads
c
c============================
c
c  If there is only one data file than type:
c
c  transform columns=!k1col,!k2col,...,!kncol,
c
c  The !kicol are the column numbers to be transformed using the settings
c  below. Commas separate each column number. The LAST integer must be
c  followed by a comma. The number of columns used must be <=48.
c
c  If there are !nfil data files then:
c
c  transform columns=!kcol,
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 or Hz, OR sr is a sampling interval with the unit determined by
c   the unit determined by the characters placed after unit=
c   If there is an 'a' in the format, then a time mark file is read after
c   the data is read for each time point.
c
c=======================================================================
c
c  file read start=!is  file read end=!ie
c
c  !is - index of first datum to read, !ie - index of last datum to read
c  If !ne = 0 the wholse data file is read.
c
c============================
c
c  Subtract column col1 from column col2=!k1col,!k2col,
c
c  !kicol is subtracted from !k2col row by row. The comma separated the
c  column numbers and the 2nd number must be followed by a comma.
c
c============================
c
c  If there is only one file:
c  The loop on the columns selected starts for this control file.
c  The data transformation commands are then read.
c  Repeat the transformation lines for each column.
c
c  If Output transformed columns in one file output=yes is  typed then
c  the output is written to a multi-column file called Rearranged.data.
c
c  If there are several files where one column from each file is selected
c  the output is written to a multi-column file called Rearranged.data.
c
c  The loop on the file columns selected starts at this point in the cnl file.
c  The data transformation commands are then read.
c  Repeat the transformation lines for each file.
c
c  Each of the list of command lines are followed by the line:
c  # */ End of column transformation switches
c
c============================
c
c  Type one of the following character strings to control the transformations
c  of the data. If neither of these two character strings are not typed in the
c  cnl file then the data will NOT be transformed.
c
c  IMPORTANT: The character strings must be in LOWER CASE
c
c=========================
c
c  column=!kr
c
c  !kr is the number of the column to be transformed in the loop on column
c
c=========================
c
c  lower level for clipping=!th
c
c  If this command is typed in the control file then a single x(t)<th is set
c  at th. If there is a consecutive string of x(t)'s<th then the string is
c  set at the moving average of the x(t)'s >= th at the beginning of the string.
c
c=========================
c
c  log
c
c  If the character string 'log' is present in the control file then the data
c  is tranformed by taking the natural log.
c
c=========================
c
c  rates of return
c
c  If the character string 'rates' is present in the control file then the
c  data is transformed into rates of return - log(x(t))-log(x(t-1)).
c
c=========================
c
c  first difference
c
c  If the character string 'difference' is present in the control file then
c  the data is differenced - x(t)-x(t-1).
c
c=========================
c
c  block scale=!sca  start at=!ks  end at=!ke
c
c  If the character string 'scal' is present in the control file then the data
c  in the block (ks,ke) is scaled by the value of !sca.
c
c=========================
c
c  detrend linear
c
c  If the character string 'detrend linear' is present in the control file
c  then the residuals of an iterative least squares fit of a line is written
c  to a file fname-Col_kc where kc is the column number selected above.
c  The iteration is controlled by a p-Value threshold for the t values of the
c  estimated coeffients that is discussed below.
c
c=========================
c
c  detrend curve
c
c  If the character string 'detrend curve' is present in the control file
c  then the residuals of an iterative least squares fit of a 1-4 th order
c  Legendre polynomials is written to a file fname-Col_kc.
c  The iteration is controlled by a p-Value threshold for the t values of the
c  estimated coeffients that is discussed below.
c
c  The Legendre polynomials are orthogonal.
c  The Legendre model is estimated in the subroutine detrend.for that outputs
c  the po < np used as the end the iterations.
c
c=========================
c
c  If you want to prewhiten the data (or the bandpassed data if bandpass=yes)
c  with an AR(np) fit type the character string
c
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 np=0.
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  Probability threshold=!rpv
c  !rpv - probability level threshold for the t statistics if the data is
c  detrended OR if the data is fitted to an AR.
c
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 leqast squares 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  If you want to trim the data at the top quantile %ptq type the character
c  string:
c
c  quantile level to trim back large outliers=%ptq
c
c  All values larger than the 100-%ptq will be trimmed down to that quantile.
c
c  If %ptq=0.0 or if the character string is not typed in the cnl file then
c  the large outliers in the data are not trimmed.
c
c  If you want to trim the data at the bottom quantile %pbq type the character
c  string:
c
c  quantile level to trim back small outliers=%pbq
c
c  All values smaller than the %ptq will be trimmed up to that quantile.
c
c  If %pbq=0.0 or if the character string is not typed in the cnl file then
c  the small outliers in the data are not trimmed.
c
c  If you want to output the trimmed values to a file called fname-co_trim.data
c  where 'fname' is the name of the data file and 'co' is the column number,
c  then type:
c
c  output trimmed data
c
c=========================
c
c  If you want to hard clip the data type the character string:
c
c  hard clip
c
c  If this character string is present in the control file then the data will
c  transformed into 1 if x(t)>= mean or -1 if x(t)< mean.
c
c=========================
c
c  If you want to bandpass filter the data to the band set below type the string
c
c  bandpass filter
c
c  No bandpass filtering is made if these words (character strings) are not
c  typed in the lines above the delimiter line that starts with the symbol #.
c
c============================
c
c  If bandpassed filter is set for the band (!rl,!ru) set below and you want
c  the low frequency (0,!rl) component to be written to the file
c  fname-co_low.data type the string 'output low band' in a line
c
c  The file fname-co_low.data will not be opened if the words output low band
c  are not in the cnl file.
c
c============================
c
c  If you type the string 'bandpass filter' to filter the data in the frequency
c  domain then set the bandlimits as follows depending on whether frequency or
c  time units are used:
c
c  Lower frequency=!rl  upper frequency=!ru
c  in units of kHz if 'frequency+k' is used, or Hz if 'frequency+h' is used
c
c  Longest period=!rl  shortest period=!ru if a unit is used
c
c  (!rl,!ru) is the bandwidth of the output data
c
c  If !rl= 0  then the lower frequency is 0.0 or the longest period is n where
c  n = !ne-!ns is the sample size.
c
c  If !ru= 0  then the normalized upper frequency is 0.5  or the
c  normalized shortest period is lb/2.
c
c=========================
c
c  If bandpassed filter is set type the word
c
c  decimate
c
c  to obtain an equally spaced output whose sampling interval is 2/!ru is the
c  sampling unit is Hz or kHz OR 2*!ru if the unit is time. If the character
c  string 'deciminate' is not typed on a line of the cnl file then the filter
c  output has the sample time interval as the raw data.
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
c   nfl - no. of observations, sr - sampling rate in kHz if freq. is
c   set, or in multiples of the time unit set in data-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===========================================
      allocatable::x,xo,ts
      parameter(mt=3,ndy=1,ny=2009)
      real*8 x(:),xo(:,:)
      real gs,rt
      integer nums(50,2),checkn,cnt,fils,kf
      character*700 name
      character*120 fn
      character*80 buf,cnl,fid
      character*50 charead,fmt,fname
      character*25 ts(:)
      character*20 par
      character*10 su
      character*2 delim
      character mk,f1o
      logical exs
c**********************
      open(9,file='Data-error.out')
c  Inquire if file exists
      inquire(file='Rundata.cnl',exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'Rundata.cnl file does not exist'
       write(9,*)'Rundata.cnl file does not exist'
       stop
      endif
c  Open control file
      open(1,file='Rundata.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 k=1,80
        buf(k:k)=' '
       enddo
c  Read line from cnl file
       read(1,'(a80)',end=2) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line of the cnl fil
     *e'
        write(9,*)'Place comment marker * after each line in the cnl fil
     *e'
        write(6,'(a70)')buf
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
c  End of reads from Rundata.cnl file=1
      enddo
c-------
c  Data control file name
      delim(1:2)='= '
      par='ile'
      cnl=charead(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,'('' Error in read of data cnl file '',a80/)')cnl
       write(9,'('' Error in read of data cnl file '',a80/)')cnl
       stop
      endif
      inquire(file=cnl,exist=exs)
      if(.not.exs) then
       write(6,'(/)')
       write(6,*)'data cnl file does not exist ',cnl
       write(9,'(/)')
       write(9,*)'data cnl file does not exist ',cnl
       stop
      endif
c  Open data control file
      open(2,file=cnl,err=3,iostat=io3,status='old')
c-------
c  Number of data files
      par='ta files'
      fils=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(6,'('' Error in read of the number of data files '',i2/)')f
     *ils
       write(9,'('' Error in read of the number of data files '',i2/)')f
     *ils
       stop
      endif
c-------
c  Check on data output in one file
      f1o='n'
      par='put'
      f1o=charead(par,name,delim,ia,ier,3)
      if(ier>0.and.(f1o/='y'.and.f1o/='Y')) then
       write(6,'('' Error in read of the data output structure '',a1/)')
     *f1o
       write(9,'('' Error in read of the data output structure '',a1/)')
     *f1o
       stop
      endif
      if(f1o=='Y') then
c  Convert to lower case y
       f1o='y'
      endif
c-------
      if(fils>1.or.f1o=='y') then
c  fn - data file name
       open(4,file='Rearranged.data')
      endif
c-------
      checkn=0
c  Loop on data files
      do kf=1,fils
       read(1,'(a80)',end=5) buf
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line of the cnl fil
     *e'
        write(9,*)'Place comment marker * after each line in the cnl fil
     *e'
        write(6,'(a70)')buf
        write(9,'(a70)')buf
        stop
       endif
c-------
c  Data file name from Rundata.cnl file
       par='ata name'
       fn=charead(par,buf,delim,ia,ier,3)
       if(ier>0) then
        write(6,*)'Error in read of data file name in Runddata.cnl ',fn
        write(9,*)'Error in read of data file name in Runddata.cnl ',fn
        stop
       endif
       inm=index(fn,'.')
       if(inm==0) then
        write(6,*)'Data file ',fn,' has no extension'
        write(9,*)'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(9,'(/)')
        write(9,*)'Input file ',fn(:inm+5),' does not exist'
        stop
       endif
c-------
      if(fils==1) then
c  Open output file
       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
       fname(nfn:nfn+3)='.out'
       open(3,file=fname(:nfn+3),err=6,iostat=io6)
      endif
c-------
c  Name of program
      write(3,'(7x,''Data Read & Transformation Program - Melvin Hinich'
     *'/)')
c  Open kf data file
       open(kf+6,file=fn,err=4,iostat=io4,status='old')
       call readhead(buf,nfl,nc,sr,fmt,iunit,kf+6,io)
       if(fils>1) then
        write(3,'(17x,''File Number '',i2)')kf
       endif
       write(3,'(1x,a79)')buf
       rewind(kf+6)
c  Trap unequal sample sizes for multiple data files
       if(nfl/=checkn.and.kf>1) then
        write(6,'(/2x,''The data files do not have the same sample sizes
     *''/)')
        write(9,'(/2x,''The data files do not have the same sample sizes
     *''/)')
        stop
       endif
       checkn=nfl
c-------
c  Allocate xo for 48 columns
       allocate(xo(0:nfl-1,0:47),stat=ibad)
       if(ibad/=0) then
        write(3,*)'Unable to allocate work space for xo in data'
        write(6,*)'Unable to allocate work space for xo in data'
        write(9,*)'Unable to allocate work space for xo in data'
        stop
       endif
c-------
       mk='a'
c  Parse format to see if there is an 'a' for time values
c  If so then set mk='m' to flag times
       ia=index(fmt,'a')
       if(ia>0) then
        mk='m'
       endif
c-------
       if(kf==1) then
        allocate(ts(0:nfl-1),stat=ibad)
        if(ibad/=0) then
         write(3,*)'Unable to allocate work space for date/time'
         write(6,*)'Unable to allocate work space for date/time'
         write(9,*)'Unable to allocate work space for date/time'
         stop
        endif
c-------
        call top(xo,fn,cnl,fils,nums,nfl,mk,ts,cnt,sr,fid,fmt,su,is,ie,n
     *,iunit,nc,nst,kf,3)
       endif
c-------
       allocate(x(0:n*nc-1),stat=ibad)
       if(ibad/=0) then
        write(3,*)'Unable to allocate work space for x in data'
        write(6,*)'Unable to allocate work space for x in data'
        stop
       endif
       call rc(x,xo,nums,ts,is,ie,iunit,n,cnt,nst,sr,su,mk,fmt,cnl,fn,nc
     *,fid,fils,f1o,kf,3)
c-------
       write(3,'(/2x,''First datum used '',i8,2x,''Last datum used'',i8,
     *2x,''Sample size = '',i8/)')is,ie,n
       if(iunit==1) then
        write(3,'(7x,''Sampling rate ='',f12.3,'' kHz'')')sr
         elseif(iunit==2) then
          write(3,'(7x,''Sampling rate ='',f12.3,'' Hz'')')sr
           else
            write(3,'(7x,''Sampling interval = '',g12.6,1x,a7/)')sr,su
       endif
       write(3,'(1x,40(''#'')/)')
       deallocate(x)
c  End of loop on data files
      enddo
c-------
c  Output data columns to Rearranged.data if fils>1 or f1o=y
      if(fils>1.or.f1o=='y') then
       write(4,'(a80)')buf
       nfd=fils-1
       ncd=nc-1
       if(mk=='m') then
        buf='(a25,2x,:,40(f17.6,1x))'
        if(iunit>0) then
         write(4,'(''rows='',i7,'' columns='',i2,'' sampling rate='',f7.
     *3,2x,''format='',a32,'' Unit='',a5)')nfl,fils,sr,buf,su
        endif
        if(iunit==0) then
         write(4,'(''rows='',i7,'' columns='',i2,'' sampling interval=1.
     *0'',2x,''format='',a32,'' Unit='',a5)')nfl,fils,buf,su
        endif
        if(f1o=='n') then
         do t=0,nfl-1
          write(4,buf)ts(t),(xo(t,kmo),kmo=0,nfd)
         enddo
        endif
        if(f1o=='y') then
         do t=0,nfl-1
          write(4,buf)ts(t),(xo(t,kmo),kmo=0,ncd)
         enddo
        endif
       endif
       if(mk/='m') then
        buf='(:,40(f17.6,1x))'
        if(iunit>0) then
         write(4,'(''rows='',i7,'' columns='',i2,'' sampling rate='',f7.
     *3,2x,''format='',a32,'' Unit='',a5)')nfl,fils,sr,buf,su
        endif
        if(iunit==0) then
         write(4,'(''rows='',i7,'' columns='',i2,'' sampling interval=1.
     *0'',2x,''format='',a32,'' Unit='',a5)')nfl,fils,buf,su
        endif
        if(f1o=='n') then
         do t=0,nfl-1
          write(4,buf)(xo(t,kmo),kmo=0,nfd)
         enddo
        endif
        if(f1o=='y') then
         do t=0,nfl-1
          write(4,buf)(xo(t,kmo),kmo=0,ncd)
         enddo
        endif
       endif
      endif
c  End of rearrange output writes
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      deallocate(xo,ts)
      stop
  1   write(6,*)'**** Error on opening Rundata.cnl ',io1
      write(9,*)'**** Error on opening Rundata.cnl ',io1
      stop
  2   write(6,'(7x,''ERROR''/)')
      write(6,*)'End of file in the read of the Rundata.cnl parameters'
      write(6,*)'Place the symbol # at the end of the parameter list in
     * the Rundata.cnl file'
      write(6,'(a80)')cnl
      write(9,*)'End of file in the read of the Rundata.cnl parameters'
      write(9,*)'Place the symbol # at the end of the parameter list in
     * the Rundata.cnl file'
      write(9,'(a80)')cnl
  3   write(6,*)'**** Error on opening data.cnl file ',io3,' kf = ',kf
      write(9,*)'**** Error on opening data.cnl file ',io3,' kf = ',kf
      stop
  4   write(6,*)'**** Error on opening data file ',io4,' kf = ',kf
      write(9,*)'**** Error on opening data file ',io4,' kf = ',kf
      stop
  5   write(6,*)'Error on read of data file on rundata.cnl reads'
      write(9,*)'Error on read of data file on rundata.cnl reads'
      stop
  6   write(6,*)'**** Error on opening main out file ',io6
      write(9,*)'**** Error on opening main out file ',io6
c-------
      stop
      end
c
c===========================================
c
      subroutine top(xo,fn,cnl,fils,nums,nfl,mk,ts,cnt,sr,fid,fmt,su,is,
     *ie,n,iunit,nc,nst,kf,io)
c********************************************************
c  Runs each control file in the master Rundata.cnl where kf - data file index
c  fn - data file name, cnt - no. of columns used, nums- column numbers used
c========================================================
      integer nums(50,2),cnt,fils,kf
      real*8 xo(nfl,48)
      real gs,rt
      character*700 name
      character*120 fn
      character*80 buf,cnl,fid
      character*50 charead,fmt
      character*25 ts(nfl)
      character*20 par
      character*10 su
      character*5 timec
      character*2 delim
      character mk,f1o
      logical exs
      intent(in) cnl,fn,fils,nfl,mk,kf,io
      intent(out) nums,cnt,sr,fid,fmt,su,is,ie,n,iunit,nc,nst
c**********************
      ip=0
      kt=0
      ia=0
      dowhile(kt<=17.and.ia==0)
       kt=kt+1
c  Clear buf
       do k=1,80
        buf(k:k)=' '
       enddo
c  Read line from cnl file
       read(2,'(a80)',end=1) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(6,*)'Place comment marker * after each line of the cnl fil
     *e'
        write(9,*)'Place comment marker * after each line in the cnl fil
     *e'
        write(6,'(a70)')buf
        write(9,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
c  End of reads for file=1
      enddo
c-------
      delim(1:2)='= '
c  Transform these data columns
      par='olumns'
      buf=charead(par,name,delim,ia,ier,3)
c  Columns numbers in nums
      call numsequence(buf,nums,cnt,ier,3)
      if(ier>0) then
       write(6,'(/7x,''ERROR-reading columns used '',i3,'' Check trailin
     *g comma'')')cnt
       write(9,'(/7x,''ERROR-reading columns used '',i3,'' Check trailin
     *g comma'')')cnt
       stop
      endif
      if(cnt>48) then
       write(6,'(/7x,''No. of columns used > 48'')')cnt
       write(9,'(/7x,''No. of columns used > 48'')')cnt
       stop
      endif
c-------
c  Sampling unit
      ik1=0
      ih1=0
      par='unit'
      buf=charead(par,name,delim,ia,ier,3)
      su=buf(:15)
      ik=index(buf,'Kh')
      ih=index(buf,'Hz')
      ikl=index(buf,'kh')
      ihl=index(buf,'hz')
      if(ik>0 .or. ikl>0) then
c  kHz
       iunit=1
       su='kHz'
        elseif(ih>0 .or. ihl>0) then
c  Hz
         iunit=2
          su='Hz'
           elseif(ik==0 .and. ih==0 .and. ik1==0 .and. ih1==0) then
c  Unit set in cnl file
            iunit=0
      endif
c-------
c  File read start
      par='start'
      is=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in read of start datum'
       write(io,*)par
       write(6,*)'Error in read of start datum'
       write(6,*)par
       write(9,*)'Error  in read of start datum'
       write(9,*)par
       stop
      endif
c  File read end
      par='end'
      ie=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(io,*)'Error in read of end datum'
       write(io,*)par
       write(6,*)'Error in read of end datum'
       write(6,*)par
       write(9,*)'Error in read of end datum'
       write(9,*)par
       stop
      endif
      if(is<1) then
       write(io,*)'IS in %.cnl is < 1 ',is
       write(6,*)'IS in %.cnl is < 1 ',is
       write(9,*)'IS in %.cnl is < 1 ',is
       stop
        elseif(ie/=0.and.ie<is) then
         write(io,*)'IE < IS in %.cnl ',is,ie
         write(6,*)'IE < IS in %.cnl ',is,ie
         write(9,*)'IE < IS in %.cnl ',is,ie
         stop
      endif
c  Check for column subtraction
      nst=0
      ik=index(name,'ubtrac')
      if(ik>0) then
c  Subtract col1 from col2
       par='ol2'
       buf=charead(par,name,delim,ia,ier,3)
       if(ier>0) then
        write(io,'(/2x,''Error in read of column subtraction''/)')
        write(io,'(3x,''Did you end the command line with col2?'')')
        write(6,'(/2x,''Error in read of column subtraction''/)')
        write(6,'(3x,''Did you end the command line with col2?'')')
        write(6,*)'Did you end the command line with col2?'
        write(9,'(/2x,''Error in read of column subtraction''/)')
        write(9,'(3x,''Did you end the command line with col2?'')')
        stop
       endif
       call numsequence(buf,nums(:49,2),nst,ier,io)
       if(nst/=2) then
        write(io,'(/2x,''Error in read of the two column numbers''/)')
        write(6,'(/2x,''Error in read of the two column numbers''/)')
        write(9,'(/2x,''Error in read of the two column numbers''/)')
        stop
       endif
      endif
c------
c  Read data file header
      call readhead(fid,nfl,nc,sr,fmt,iunit,kf+6,io)
c  Trap no. of columns nc
      if(nc<1) then
       write(io,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(6,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       write(9,'(/7x,''Number of columns '',i3,'' < 1'')')nc
       stop
      endif
      if(cnt>nc) then
       write(io,'(/7x,''No. of columns set in the cnl file '',i2,'' < no
     *. of columns in the data file'')')cnt,nc
       write(6,'(/7x,''No. of columns set in the cnl file '',i2,'' < no.
     * of columns in the data file'')')cnt,nc
       write(9,'(/7x,''No. of columns set in the cnl file '',i2,'' < no.
     * of columns in the data file'')')cnt,nc
       stop
      endif
      if(cnt==0) then
       cnt=nc
      endif
c  Trap file length
      if(ie>nfl) then
       write(io,*)'IE ',ie, ' > N in data file',nfl
       write(6,*)'IE ',ie, ' > N in data file',nfl
       write(9,*)'IE ',ie, ' > N in data file',nfl
       stop
      endif
c-------
      if(ie==0) then
       ie=nfl
      endif
      n=ie-is+1
      return
  1   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(9,*)'End of file in the read of the cnl parameters'
      write(9,*)'Place the symbol # at the end of the parameter list in
     * the cnl file'
      write(9,'(a80)')cnl
      stop
      end
c
c===========================================
c
      subroutine rc(x,xo,nums,ts,is,ie,iunit,n,cnt,nst,sr,su,mk,fmt,cnl,
     *fn,nc,fid,fils,f1o,kf,io)
c********************************************************
c  Read data & loop on columns for kf data file
c  cnt - no. of columns used, fn - data file name, nums- column numbers used
c========================================================
      allocatable::ao,c,e,ik,r,rt,s,w,y,z
      real*8 x(n,nc),xo(n,48),y(:),z(:),r(:),c(:),e(:,:),s(:),w(:),del,a
     *m,sd,c3,c4,max,min
      real fl,fu
      complex ao(:),rt(:)
      integer ik(:),nums(50,2),cnt,nst,p,fils,kf,t
      character*700 name
      character*120 fn
      character*80 buf,cnl,fid
      character*50 fmt,fname,fout
      character*25 ts(n)
      character*20 par
      character*10 su
      character*3 cntr
      character*2 co,dc,delim
      character bpf,f1o,mk,tr
      intent(in) nums,is,ie,iunit,n,cnt,nst,sr,su,mk,fid,fmt,cnl,fn,nc,f
     *ils,f1o,kf,io
      intent(out) xo,ts
c**********************
c  Read data
      ic=index(fmt,'*')
      if(ic>0) then
c  * read format
       call dats(x,is,ie,n,nc,kf+6,io)
      endif
      if(ic==0) then
c  Formatted data read with or without a character time/date column
       call datc(x,ts,fmt,mk,is,ie,n,nc,kf+6,io)
      endif
      if(nst==2) then
       nc1=10*nums(49,1)+nums(49,2)
       nc2=10*nums(50,1)+nums(50,2)
c  Subtract column nc1 from nc2 =>nc2
       do t=1,n
        x(t,nc2)=x(t,nc2)-x(t,nc1)
       enddo
      endif
c------
c  Allocate space
      ia=n-1
      allocate(e(0:ia,0:1),r(0:ia),s(0:ia),w(0:ia),y(0:ia),z(0:ia),ik(0:
     *ia),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays in rc'
       write(io,'(/2x,''Total array size called = '',i12/)')7*n
       write(6,*)'Unable to allocate work space for arrays in rc'
       write(6,'(/2x,''Total array size called = '',i12/)')7*n
       write(9,*)'Unable to allocate work space for arrays in rc'
       write(9,'(/2x,''Total array size called = '',i12/)')7*n
       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
c-------
      fname=fn(ib:ia)
      fout=fn(ib:ia)
      nfn=ia-ib+1
      delim(1:2)='= '
c  Loop on columns
      write(io,'(7x,''Data Columns in the File = '',i2)')nc
      do kr=1,cnt
c-------
c  Clear name
       do t=1,700
        name(t:t)=' '
       enddo
c  Column no in the loop
       kcol=10.0*nums(kr,1)+nums(kr,2)
       ip=0
       kt=0
       ia=0
       dowhile(kt<=17.and.ia==0)
        kt=kt+1
c  Clear buf
        do t=1,80
         buf(t:t)=' '
        enddo
c  Read transformation lines
        read(2,'(a80)',end=4) buf
        ia=index(buf,'#')
        ig=index(buf,'*')
        if(ig==0) then
         write(6,*)'Place comment marker * after each line of the cnl fi
     *le kf= ',kf
         write(6,'(2x,''Column = '',i2/)')kr
         write(9,*)'Place comment marker * after each line in the cnl fi
     *le kf= ',kf
         write(9,'(2x,''Column = '',i2/)')kr
c-------
         write(6,'(a70)')buf
         write(9,'(a70)')buf
         stop
        endif
        ig=ig-1
        name(ip+1:ip+ig)=buf(:ig)
        ip=ip+ig
c  End of reads of transformataion lines
       enddo
c-------
c  Column number
       nc1=10*nums(kr,1)+nums(kr,2)
       par='umber'
       knum=numread(par,name,delim,ia,ier,io)
       if(knum/=nc1) then
        write(6,'(2x,''Column number in cnl '',i2,'' is not equal to the
     *column index '',i2/)')knum,nc1
        write(9,'(2x,''Column number in cnl '',i2,'' is not equal to the
     *column index '',i2/)')knum,nc1
        stop
       endif
c-------
c  Bandpass filter
       bpf='n'
       gl=0.0
       gu=0.5
       ia=index(name,'andpa')
       if(ia>0) then
        bpf='f'
       endif
       dc='nd'
       ia=index(name,'cimate')
       if(ia>0) then
c  Decimate switch
        dc='de'
       endif
c-------
       if(iunit>0.and.bpf=='f') then
c  Lower frequency
        par='wer freq'
        fl=rnumread(par,name,delim,ia,ier,io)
        if(ier==1) then
         write(io,'(/'' Type lower freq rather than longest period in se
     *tting the band in %.cnl'')')
         write(6,'(/'' Type lower freq rather than longest period in set
     *ting the band in %.cnl'')')
         write(9,'(/'' Type lower freq rather than longest period in set
     *ting the band in %.cnl'')')
         return
        endif
c  Upper frequency
        par='per freq'
        fu=rnumread(par,name,delim,ia,ier,io)
        if(ier==1) then
         write(io,'(/'' Type upper freq rather than shortest period in s
     *etting the band in %.cnl'')')
         write(6,'(/'' Type upper freq rather than shortest period in se
     *tting the band in %.cnl'')')
         write(9,'(/'' Type upper freq rather than shortest period in se
     *tting the band in %.cnl'')')
         return
        endif
c-------
        if(ier>0) then
         write(io,'(/7x,''Set band limits for bandpass filtering''/)')
         write(io,'(7x,''Lower & upper bandlimits = '',2(f12.3,2x))')fl,
     *fu
         write(6,'(/7x,''Set band limits for bandpass filtering''/)')
         write(6,'(7x,''Lower & upper bandlimits = '',2(f12.3,2x))')fl,f
     *u
         write(9,'(/7x,''Set band limits for bandpass filtering''/)')
         write(9,'(7x,''Lower & upper bandlimits = '',2(f12.3,2x))')fl,f
     *u
         return
        endif
c-------
        if(fu>0.0.and.fl>fu) then
         write(io,*)'Upper bandlimit in %.cnl ',fu,' < lower bandlimit '
     *,fl
         write(6,*)'Upper bandlimit in %.cnl ',fu,' < lower bandlimit ',
     *fl
         write(9,*)'Upper bandlimit in %.cnl ',fu,' < lower bandlimit ',
     *fl
         if(iunit==1) then
          write(io,*)'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(io,*) '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'
         endif
         return
        endif
c-------
c  Convert to normalized values (0,.5)
        gl=fl/sr
        if(fu>0.0) then
         gu=fu/sr
          else
           gu=0.5
        endif
       endif
c-------
       if(iunit==0.and.bpf=='f') then
c  Longest period
        par='gest period'
        fl=rnumread(par,name,delim,ia,ier,io)
        if(ier==1) then
         write(io,'(/'' 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 set
     *ting the band in the cnl file'')')
         write(9,'(/'' Type longest period rather than lower freq in set
     *ting the band in the cnl file'')')
         return
        endif
c  Shortest period
        par='test period'
        fu=rnumread(par,name,delim,ia,ier,io)
        if(ier==1) then
         write(io,'(/'' Type shortest period rather than upper freq in s
     *etting the band in %.cnl'')')
         write(6,'(/'' Type shortest period rather than upper freq in se
     *tting the band in %.cnl'')')
         write(9,'(/'' Type shortest period rather than upper freq in se
     *tting the band in %.cnl'')')
         return
        endif
        if(fl>0.0.and.fu>0.0.and.fl<fu) then
         write(io,*)'Longest period in %.cnl ',fl,' < shortest period ',
     *fu
         write(6,*)'Longest period in %.cnl ',fl,' < shortest period ',f
     *u
         write(9,*)'Longest period in %.cnl ',fl,' < shortest period ',f
     *u
         return
        endif
c-------
c  Trap fl=0 for low band period
        if(fl<=0.0) then
         gl=0.0
          elseif(fl>0.0) then
c  Invert period to frequency
           gl=1.0/fl
        endif
c  Convert to normalized values (0,.5)
        gl=gl*sr
        if(fu>0.0) then
         gu=(1.0/fu)*sr
          else
           gu=0.5
        endif
       endif
c-------
c  Bandpass
       if(bpf=='f') then
c  Trap top band
        if(gu>0.5) then
         write(io,'(/2x,''Normalized upper bandlimit = '',f10.4,''  > 0.
     *5   Column '',i3/)')gu,kcol
         write(6,'(/2x,''Normalized upper bandlimit = '',f10.4,''  > 0.5
     *   Column '',i3/)')gu,kcol
         write(9,'(/2x,''Normalized upper bandlimit = '',f10.4,''  > 0.5
     *   Column '',i3/)')gu,kcol
         if(iunit==1) then
          write(io,*) 'Check to see if the (fl,fu) units are kHz'
          write(6,*) 'Check to see if the (fl,fu) units are kHz'
          write(9,*) 'Check to see if the (fl,fu) units are kHz'
           elseif(iunit==2) then
            write(io,*) 'Check to see if the (fl,fu) units are Hz'
            write(6,*) 'Check to see if the (fl,fu) units are Hz'
            write(9,*) 'Check to see if the (fl,fu) units are Hz'
             else
              write(io,*) 'Check to see if the (fl,fu) units are ',su
              write(6,*) 'Check to see if the (fl,fu) units are ',su
              write(9,*) 'Check to see if the (fl,fu) units are ',su
         endif
         write(io,*) ' Sampling unit ',su,' Sampling rate/interval ',sr
         write(io,*) '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
         return
        endif
       endif
c-------
       tr='n'
       ia=index(name,'utput')
       if(ia>0) then
c  Output low frequency components
        tr='t'
       endif
c-------
       if(nc>1.and.fils==1) then
        write(io,'(/12x,27(''=''))')
        write(io,'(17x,''| Column = '',i2,'' |'')')kcol
        write(io,'(12x,27(''=''))')
       endif
c-------
c  Column kr name
       if(kcol<10.and.fils==1.and.f1o=='n') then
        co(1:1)=achar(kcol+48)
        co(2:2)=' '
        fname(nfn:nfn+4)='-Col_'
        fname(nfn+5:nfn+5)=co
        fname(nfn+6:nfn+10)='.data'
        fout(nfn:nfn)='-'
        fout(nfn+1:nfn+2)=co
       endif
       if(kcol>9.and.fils==1.and.f1o=='n') then
        co(1:1)='1'
        co(2:2)=achar(kcol+38)
        fname(nfn:nfn+4)='-Col_'
        fname(nfn+5:nfn+6)=co
        fout(nfn:nfn)='-'
        fout(nfn+1:nfn+2)=co
       endif
       if(fils==1.and.f1o=='n') then
c  Open fname-co_#.data
        close(47+kr)
        open(47+kr,file=fname(:nfn+11),err=1,iostat=io1)
        write(47+kr,'(''Column '',a2,2x,a77)')co,fid
        if(bpf=='f'.and.tr=='t') then
c  Open fname-co_low.data
         fout(nfn+3:nfn+11)='_low.data'
         close(97+kr)
         open(97+kr,file=fout(:nfn+11),err=2,iostat=io2)
        endif
       endif
c  Shift kcol column to y
       do t=1,n
        y(t-1)=x(t,kcol)
       enddo
c-------
       cntr='nnnn'
       pcu=0.0
       pcl=0.0
       ia=index(name,'ates')
       if(ia>0) then
c  Rates of return
        cntr(1:1)='r'
       endif
c-------
       cth=-1.e-17
       ia=index(name,'clip')
       if(ia>0) then
c  Level for clipping
        par='ping'
        cth=rnumread(par,name,delim,ia,ier,io)
        if(ier>0) then
         write(io,'(/2x,''ERROR reading the clip threshold ''/)')cth
         write(6,'(/2x,''ERROR reading the clip threshold ''/)')cth
         write(9,'(/2x,''ERROR reading the clip threshold ''/)')cth
         stop
        endif
c  Data clipped at level cth
        write(io,'(/7x,''Data is Clipped at the Lower Level = '',f10.3)'
     *)cth
c  Statistics of the data
        call dstat(n,y,am,sd,c3,c4,max,min,io)
        write(io,'(60(''*''))')
        write(io,'(2x,''Statistics of the Data before Clipping'')')
        call dwrstat(io,am,sd,c3,c4,max,min)
        call dataclp(y,dble(cth),n,io)
       endif
c-------
       nsc=0
       ia=index(name,'cale')
       if(ia>0) then
c  Scale for block scaling
        nsc=1
        par='cale'
        sca=rnumread(par,name,delim,ia,ier,io)
        if(ier>0) then
         write(6,'(/2x,''ERROR reading the frame scale value ''/)')
         write(9,'(/2x,''ERROR reading the frame scale value ''/)')
         stop
        endif
       endif
       if(nsc==1) then
        par='t at'
        ks=numread(par,name,delim,ia,ier,io)
        if(ier>0) then
         write(6,'(/2x,''ERROR reading 1st datum to be scaled ''/)')
         write(9,'(/2x,''ERROR reading 1st datum to be scaled ''/)')
         stop
        endif
        par='d at'
        ke=numread(par,name,delim,ia,ier,io)
        if(ier>0) then
         write(6,'(/2x,''ERROR reading 1ast datum to be scaled ''/)')
         write(9,'(/2x,''ERROR reading 1ast datum to be scaled ''/)')
         stop
        endif
        if(ke==0) then
         ke=n
        endif
c-------
c Trap ks & ke
        if(ks>ke) then
         write(6,'(2x,''ERROR - start '',i7,'' > end '',i8/)')ks,ke
         write(9,'(2x,''ERROR - start '',i7,'' > end '',i8/)')ks,ke
        endif
        if((ke-ks+1)>n) then
         write(6,'(2x,''ERROR: end - start +1 '',i7,'' > n '',i8/)')ke-k
     *s+1,n
         write(9,'(2x,''ERROR: end - start +1 '',i7,'' > n '',i8/)')ke-k
     *s+1,n
         stop
        endif
c-------
        write(io,'(/2x,''Data Block '',i7,'' - '',i8,'' is scaled by '',
     *f10.5)')ks,ke,sca
        do t=ks-1,ke-1
         y(t)=sca*y(t)
        enddo
       endif
c-------
       ia=index(name,'og')
       if(ia>0) then
c  Log data
        cntr(1:1)='g'
       endif
       ia=index(name,'erence')
       if(ia>0) then
c  First differences
        cntr(1:1)='f'
       endif
       ia=index(name,'ard cli')
       if(ia>0) then
c  Hard clip around mean
        cntr(3:3)='c'
       endif
       ia=index(name,'end linear')
       if(ia>0) then
c  Detrend using a least squares fit of centered line
        cntr(3:3)='l'
       endif
       ia=index(name,'end curve')
       if(ia>0) then
c  Detrend using a least squares fit of centered quadratic polynomial
        cntr(3:3)='d'
       endif
c-------
c  Prewhiten by an AR(p) fit
       p=0
       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(ib:),delim,i1,ier,io)
          elseif(ic>0) then
           par='r'
           delim(1:2)='()'
           p=numread(par,name(ic:),delim,i1,ier,io)
        endif
       endif
       if(ia>0.and.ier>0) then
        write(6,'(/2x,''ERROR reading the AR lag p for prewhitening '',i
     *3/)')p
        write(9,'(/2x,''ERROR reading the AR lag p for prewhitening '',i
     *3/)')p
        return
       endif
c  Trap p
       if(p<0) then
        write(io,*)'AR parameter ',p,' < 0. p set = 0'
        write(6,*)'AR parameter ',p,' < 0. p set = 0'
        write(9,*)'AR parameter ',p,' < 0. p set = 0'
        p=0
         elseif(p>(n/4)) then
          write(io,*)'AR parameter ',p,' > n/4. p set = 0'
          write(6,*)'AR parameter ',p,' > n/4. p set = 0'
          write(9,*)'AR parameter ',p,' > n/4. p set = 0'
          p=0
       endif
       delim(1:2)='= '
       pv=1.0
       if(p>0.or.cntr(3:3)=='l'.or.cntr(3:3)=='d') then
c  Probability threshold for AR or Legendre t values
        par='ty threshold'
        pv=rnumread(par,name,delim,ia,ier,io)
c  Trap pv
        if(pv<=0.0) then
         write(io,*)'Prob threshold pv ',pv,' < 0'
         write(io,*)'Set pv = 0 if you do not want to prewhiten'
         write(6,*)'Prob threshold pv ',pv,' < 0'
         write(6,*)'Set pv = 0 if you do not want to prewhiten'
         write(9,*)'Prob threshold pv ',pv,' < 0'
         write(9,*)'Set pv = 0 if you do not want to prewhiten'
         stop
        endif
        if(pv>1) then
         write(io,*)'Prob threshold pv ',pv,' > 1'
         write(6,*)'Prob threshold pv ',pv,' > 1'
         stop
        endif
       endif
c-------
       if(p>0) then
        allocate(ao(0:p),c(0:p-1),rt(0:p-1),stat=ibad)
       endif
       if(p==0) then
        allocate(ao(0:0),c(0:0),rt(0:0),stat=ibad)
       endif
       if(ibad/=0) then
        write(io,*)'Unable to allocate work space for AR arrays in rc'
        write(io,'(/2x,''Total array size called = '',i12/)')3*p
        write(6,*)'Unable to allocate work space for AR arrays in rc'
        write(6,'(/2x,''Total array size called = '',i12/)')3*p
        write(9,*)'Unable to allocate work space for AR arrays in rc'
        write(9,'(/2x,''Total array size called = '',i12/)')3*p
        return
       endif
c-------
c  Trap bandpass
       if((cntr(3:3)=='l'.or.cntr(3:3)=='d').and.bpf=='f') then
        write(io,'(/2x,''You can not bandpass & detrend''/)')
        write(6,'(/2x,''You can not bandpass & detrend''/)')
        write(9,'(/2x,''You can not bandpass & detrend''/)')
        stop
       endif
       if(cntr(3:3)=='c'.and.bpf=='f') then
        write(io,'(/2x,''You can not bandpass & hard clip''/)')
        write(6,'(/2x,''You can not bandpass & hard clip''/)')
        write(9,'(/2x,''You can not bandpass & hard clip''/)')
        stop
       endif
c  Trap hard clipping
       if(cntr(1:1)=='c'.and.(cntr(3:3)=='l'.or.cntr(3:3)=='d')) then
        write(io,'(/2x,''You can not hard clip & detrend''/)')
        write(6,'(/2x,''You can not hard clip & detrend''/)')
        write(9,'(/2x,''You can not hard clip & detrend''/)')
        stop
       endif
       if(cntr(1:1)=='c'.and.p>0) then
        write(io,'(/2x,''You can not hard clip & fit an AR model''/)')
        stop
       endif
c-------
       ia=index(name,'rim')
       if(ia>0) then
c  Trim tailes
        cntr(2:2)='t'
c  Trap hard clipping
        if(cntr(1:1)=='c'.and.cntr(2:2)=='t') then
         write(io,'(/2x,''You can not hard clip & trim''/)')
         write(6,'(/2x,''You can not hard clip & trim''/)')
         write(9,'(/2x,''You can not hard clip & trim''/)')
         stop
        endif
c  Quantile for large outliers
        par='rge outliers'
        pcu=rnumread(par,name,delim,ip,ier,io)
        if(ier>0) then
         write(io,'(2x,''Error in read of large outliers quantile''/)')
         write(6,'(2x,''Error in read of large outliers quantile''/)')
         write(9,'(2x,''Error in read of large outliers quantile''/)')
         return
        endif
c  Quantile for small outliers
        par='all outliers'
        if(ier>0) then
         write(io,'(2x,''Error in read of small outliers quantile''/)')
         write(6,'(2x,''Error in read of small outliers quantile''/)')
         write(9,'(2x,''Error in read of small outliers quantile''/)')
         return
        endif
        pcl=rnumread(par,name,delim,ip,ier,io)
       endif
       ia=index(name,'immed data')
c  Output trimmed data file
       if(ia>0.and.cntr(2:2)=='t') then
        cntr(2:2)='o'
c  Open fname-co#_trim.data
        fout(nfn+2:nfn+15)='_trimmed.data'
        close(137)
        open(137+kr,file=fout(:nfn+15),err=3,iostat=io3)
       endif
       call rn(y,r,s,w,xo,z,ts,ao,c,e,rt,ik,gl,gu,sr,su,n,p,bpf,cntr,fmt
     *,mk,tr,dc,pcu,pcl,pv,f1o,kr,nc,iunit,fils,kf,io)
       call wrh(su,n,p,pcu,pcl,pv,th,bpf,fl,fu,iunit,io)
       deallocate(ao,c,rt)
c  End of column loop
      enddo
c-------
      deallocate(e,ik,r,s,w,y,z)
      return
  1   write(6,*)'**** Error opening data output',io1
      write(9,*)'**** Error opening data output',io1
      stop
  2   write(6,*)'**** Error opening low frequency output',io2
      write(9,*)'**** Error opening low frequency output',io2
      stop
  3   write(6,*)'**** Error opening trim output',io3
      write(9,*)'**** Error opening trim output',io3
      stop
  4   write(6,'(7x,''ERROR in transformation reads- Col = '',i2/)')kr
      write(6,*)'End of file in the read of the data cnl parameters'
      write(6,*)'Place the symbol # at the end of the parameter list in
     *the data cnl file kf= ',kf
      write(6,'(a80)')cnl
      write(9,*)'End of file in the read of the data cnl parameters'
      write(9,*)'Place the symbol # at the end of the parameter list in
     *the data cnl file kf= ',kf
      write(9,'(a80)')cnl
      stop
      end
c
c===========================================
c
      subroutine rn(y,r,s,w,xo,z,ts,ao,c,e,rt,ik,fl,fu,sr,su,n,p,bpf,cnt
     *r,fmt,mk,tr,dc,pcu,pcl,pv,f1o,kr,nc,iunit,fils,kf,io)
c********************************************************
c  Compute data tranformations for kr column & kf data, n=nu in call
c  Input: y - data, su - unit control, sr - sampling rate, fl - lower
c  freqency/period, fu - upper frequency/period, kf - data file
c  cntr(1:1)=r -rates, g -log, f -1st difference, c(2:2)=t -trim
c  cntr(3:3)=c -hard clip, l -linear detrend, d -orthogonal polynomial detrend
c========================================================
      allocatable::to
      real*8 y(n),r(n),c(p),e(n,2),s(n),w(n),xo(n,48),z(n),a(4),co(4),vf
     *rac(1),am,sd,c3,c4,center,max,min,ren,tn,wl,wu
      real frac(1)
      integer ik(n),fils,kf,kr,n,nc,iunit,p,po,q,t,tru,trd
      complex ao(p+1),rt(p)
      character*50 fmt
      character*25 ts(n),to(:)
      character*10 su
      character*3 cntr
      character*2 dc
      character bpf,f1o,mk,tr
      intent(in) y,ts,n,fl,fu,sr,bpf,cntr,mk,tr,fmt,su,p,pcu,pcl,pv,kr,n
     *c,dc,f1o,iunit,fils,kf,io
      intent(out) xo
c**********************
      allocate(to(0:0),stat=ibad)
      if(ibad/=0) then
       write(3,*)'Unable to allocate work space for xo in data'
       write(6,*)'Unable to allocate work space for xo in data'
       write(9,*)'Unable to allocate work space for xo in data'
       stop
      endif
c-------
c  Compute statistics for input data
      call dstat(n,y,am,sd,c3,c4,max,min,io)
      write(io,'(/60(''*''))')
      write(io,'(7x,''Descriptive Statistics of Data - N = '',i9)')n
      call dwrstat(io,am,sd,c3,c4,max,min)
c-------
      if(cntr(1:1)=='r') then
       write(io,'(/7x,''Data Transformed into Rates of Return'')')
      endif
      if(cntr(1:1)=='g') then
       write(io,'(/7x,''Data is Logged using the Natural log'')')
      endif
      if(cntr(1:1)=='f') then
       write(io,'(/7x,''Data is First Differenced'')')
      endif
      if(cntr(3:3)=='c') then
       write(io,'(/7x,''Data is Clipped to 1/-1 Values'')')
      endif
      if(bpf=='f'.and.cntr(3:3)/='d') then
c  Bandpass filter
       write(io,'(/1x,''Data Filtered using an FFT to Remove the Low & H
     *igh Frequencies Outside the Passband'')')
      endif
c-------
      if(cntr(1:1)=='r') then
c  Rate of return y =>r
       call returns(y,r,n,io)
c  Statistics of returns
       call dstat(n,r,am,sd,c3,c4,max,min,io)
       write(io,'(60(''*''))')
       write(io,'(2x,''Descriptive Statistics of the Returns'')')
       call dwrstat(io,am,sd,c3,c4,max,min)
      endif
      if(cntr(1:1)=='f') then
c  Difference y =>r
       call difference(y,r,n,io)
c  Statistics of differences
       call dstat(n,r,am,sd,c3,c4,max,min,io)
       write(io,'(60(''*''))')
       write(io,'(2x,''Descriptive Statistics of the First Differences''
     *)')
       call dwrstat(io,am,sd,c3,c4,max,min)
      endif
      if(cntr(1:1)=='g') then
c  Log y =>r
       do t=1,n
        if(y(t)<1.e-17) then
         write(io,'(/7x,''Datum '',f12.5,'' at time '',i6,'' < 1.e-17'')
     *')y(t),t
         write(6,'(/7x,''Datum '',f12.5,'' at time '',i6,'' < 1.e-17'')'
     *)y(t),t
         return
        endif
        r(t)=dlog(y(t))
       enddo
c  Statistics of logs
       call dstat(n,r,am,sd,c3,c4,max,min,io)
       write(io,'(60(''*''))')
       write(io,'(2x,''Descriptive Statistics of the Logs'')')
       call dwrstat(io,am,sd,c3,c4,max,min)
      endif
c-------
      if((cntr(1:1)=='n').and.(cntr(2:2)/='t'.or.cntr(2:2)/='o')) then
c  y =>r for untransformed or untrimmed y scaled or not
       do t=1,n
        r(t)=y(t)
       enddo
      endif
c-------
      po=0
      if(cntr(3:3)=='d') then
c  Detrend r with orthogonal Legendre polynomials
       call ddetrend(r,w,co,ik,n,pv,po,io)
c  Residuals w =>r
       do t=1,n
        r(t)=w(t)
       enddo
      endif
c-------
      if(cntr(3:3)=='l') then
c  Linear trend in array e
       ren=dfloat(n)
       center=0.5d0-1.d0/(2.d0*ren)
       do t=1,n
        tc=dfloat(t+1)/ren
        e(t,1)=tc-center
        e(t,2)=0.d0
       enddo
c  Least squares fit of r with a linear trend =>w
       call olsrun(r,e,w,n,1,ier,io)
       if(ier>0) then
        write(io,*)'Error in olsrun ',ier
        write(6,*)'Error in olsrun ',ier
        write(9,*)'Error in olsrun ',ier
        stop
       endif
c  Residuals w =>r
       do t=1,n
        r(t)=w(t)
       enddo
      endif
c-------
      if(p>0) then
c  Least squares AR(p) fit of r - residuals =>w, sigmas =>z
       call dfar(r,c,z,w,ik,r2,sig,n,p,po,pv,ier,io)
       if(ier>0) then
        write(io,*)'Error in subroutine dar ',ier
        write(6,*)'Error in subroutine dfar ',ier
        stop
       endif
       if(po>0) then
c  Statistics of residuals
        call dstat(n,w,am,sd,c3,c4,max,min,io)
c-------
c  White test with exponent exp
        exp=0.33
        call dport(w,n,exp,ct,io)
        write(io,'(60(''=''))')
        write(io,'(7x,''Descriptive Statistics of Residuals'')')
        call dwrstat(io,am,sd,c3,c4,max,min)
        write(io,'(/7x,''Whiteness test p-value = '',f7.3/)')ct
        write(io,'(60(''*''))')
        call outar(c,z,ao,rt,ik,r2,sig,su,sr,pv,p,po,iunit,io)
c  Residuals w =>r
        do t=1,n
         r(t)=w(t)
        enddo
       endif
      endif
c-------
      if(cntr(2:2)=='t'.or.cntr(2:2)=='o') then
c  Trim r =>s if p=0 or w =>r =>s if p>0
       deallocate(to)
       allocate(to(0:n-1),stat=ibad)
       if(ibad/=0) then
        write(3,*)'Unable to allocate work space for to in data'
        write(6,*)'Unable to allocate work space for to in data'
        write(9,*)'Unable to allocate work space for to in data'
        stop
       endif
       if(p==0) then
        call dtrim(r,s,z,ts,to,n,pcu,pcl,wu,wl,tru,trd,io)
       endif
       if(p>0) then
        call dtrim(w,s,z,ts,to,n,pcu,pcl,wu,wl,tru,trd,io)
       endif
c-------
c  Statistics of the trimmed data
       call dstat(n,s,am,sd,c3,c4,max,min,io)
       write(io,'(2x,''Descriptive Statistics of Trimmed Data'')')
       call dwrstat(io,am,sd,c3,c4,max,min)
       write(io,'(/2x,''Lower & upper trim quantiles = '',2(f12.3,2x),''
     * No. Trimmed from Top & Bottom = '',2(i8,2x))')wl,wu,tru,trd
      endif
c-------
      if(cntr(2:2)=='o') then
c  Sort trimmed data in z for output of trimmed
       ntr=tru+trd
       call dsort(ik,z,ntr,1,frac,vfrac,'S',ier)
       if(ier>0) then
        write(io,'(/7x,''Error in Sort in Rn''/)')
        write(6,'(/7x,''Error in Sort in Rn''/)')
        stop
       endif
       write(137+kr,'(i7,'' Sorted Outlier Data'')')ntr
       if(mk=='m') then
        do t=1,ntr
         write(137+kr,'(f17.6,2x,a25)')z(ik(t)),to(ik(t)-1)
        enddo
       endif
       if(mk/='m') then
        do t=1,tru+trd
         write(137+kr,'(i7,1x,f17.6)')ik(t),z(ik(t))
        enddo
       endif
      endif
c-------
      if(bpf=='f'.and.cntr(3:3)/='l'.and.cntr(3:3)/='d') then
c  Bandpass filter if not detrended
       call dbpfilter(r,w,z,fl,fu,n,tr,io)
       nn=n
       ndec=0
       if(dc=='de') then
c  Decimate z =>r
        if(iunit>0) then
         ndec=nint(1.0/(2.0*fu*sr))
          elseif(iunit==0) then
           ndec=nint(1.0/(2.0*fu/sr))
        endif
        nn=n/ndec
        do t=1,nn
         r(t)=z(t*ndec)                                                                              i
        enddo
       endif
c-------
c  Compute statistics of transformed & bandpassed data
       call dstat(nn,r,am,sd,c3,c4,max,min,io)
       write(io,'(60(''*''))')
c-------
       write(io,'(2x,''Descriptive Statistics of Transformed & Bandpasse
     *d Sampled Data'')')
       call dwrstat(io,am,sd,c3,c4,max,min)
       if(tr=='t'.and.fils==1) then
c  Output low band
        if(fl>0.0) then
         write(97+kr,'('' Filtered Low Band Signal'')')
        endif
        if(mk=='m'.and.fils==1) then
         if(iunit>0) then
          write(97+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate='
     *',f7.3,2x,''format=(f17.6,2x,a25) Unit='',a5)')n,sr,su
         endif
         if(iunit==0) then
          write(97+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interv
     *al='',f7.3,2x,''format=(f17.6,2x,a25) Unit='',a5)')n,sr,su
         endif
         do t=1,n
          write(97+kr,'(f17.6,2x,a25)')z(t),ts(t)
         enddo
        endif
        if(mk/='m'.and.fils==1) then
         if(iunit>0) then
          write(97+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate='
     *',f7.3,2x,''format=(f17.6) Unit='',a5)')n,sr,su
         endif
         if(iunit==0) then
          write(97+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interv
     *al='',f7.3,2x,''format=(f17.6) Unit='',a5)')n,sr,su
         endif
         do t=1,n
          write(97+kr,'(f17.6)')z(t)
         enddo
        endif
       endif
       return
      endif
c  End of bandpass filtering
c-------
      if(p>0) then
c  AR fit
       if(mk/='m'.and.fils==1.and.f1o=='n') then
        if(iunit>0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate=''
     *,f7.3,2x,''format=(f17.9)  Unit='',a5)')n,sr,su
        endif
        if(iunit==0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interva
     *l='',f7.3,2x,''format=(f17.9)  Unit='',a5)')n,sr,su
        endif
        do t=1,n
          write(47+kr,'(f17.9)')w(t)
        enddo
       endif
       if(mk=='m'.and.fils==1.and.f1o=='n') then
        if(iunit>0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate=''
     *,f7.3,2x,''format=(f17.9,2x,a25) Unit='',a5)')n,sr,su
        endif
        if(iunit==0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interva
     *l='',f7.3,2x,''format=(f17.9,2x,a25) Unit='',a5)')n,sr,su
        endif
        do t=1,n
         write(47+kr,'(f17.9,2x,a25)')w(t),ts(t)
        enddo
       endif
c  w=>xo(t,kf)
       do t=1,n
        xo(t,kf)=w(t)
       enddo
      endif
c-------
      if(cntr(3:3)=='c') then
c  Hardclip y =>r
       call hardclip(y,r,n,dble(am),io)
      endif
c-------
      if(p==0) then
c  No AR fit
       if(mk/='m'.and.fils==1.and.f1o=='n') then
        if(iunit>0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate=''
     *,f7.3,2x,''format=(f17.9) Unit='',a5)')n,sr,su
        endif
        if(iunit==0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interva
     *l='',f7.3,2x,''format=(f17.9) Unit='',a5)')n,sr,su
        endif
c-------
        if(cntr(2:2)/='t') then
         do t=1,n
          write(47+kr,'(f17.9)')r(t)
         enddo
        endif
        if(cntr(2:2)=='t'.or.cntr(2:2)=='o') then
c  Output trimmed data
         do t=1,n
          write(47+kr,'(f17.9)')s(t)
         enddo
        endif
       endif
c-------
       if(mk=='m'.and.fils==1.and.f1o=='n') then
        if(iunit>0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling rate=''
     *,f7.3,2x,''format=(f17.9,2x,a25) Unit='',a5)')n,sr,su
        endif
        if(iunit==0) then
         write(47+kr,'(''rows='',i7,'' columns=1'',2x,''sampling interva
     *l='',f7.3,2x,''format=(f17.9,2x,a25) Unit='',a5)')n,sr,su
        endif
        if(cntr(2:2)/='t'.and.cntr(2:2)/='o') then
c  No trimming
         do t=1,n
          write(47+kr,'(f17.9,2x,a25)')r(t),ts(t)
         enddo
        endif
c  Output trimmed data
        if(cntr(2:2)=='t'.or.cntr(2:2)=='o') then
         do t=1,n
          write(47+kr,'(f17.9,2x,a25)')s(t),ts(t)
         enddo
        endif
       endif
c  r=>xo(t,kf) if f1o='n' or r=>xo(t,kr) if f1o='y'
       if(f1o=='n') then
        do t=1,n
         xo(t,kf)=r(t)
        enddo
       endif
       if(f1o=='y') then
        do t=1,n
         xo(t,kr)=r(t)
        enddo
       endif
      endif
c-------
      deallocate(to)
      return
      end
c
c===========================================
c
      subroutine datc(x,ts,fmt,mk,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*8 x(n,nc)
      character*50 fmt
      character*25 ts(n)
      character mk
      intent(in) fmt,mk,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(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,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(9,*)'**** eof in read from data file '
      write(9,*)'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(9,*) '** Error in read of data file ',ko
      write(9,*) 'Check format',fmt
      stop
      end
c
c===========================================
c
      subroutine dats(x,is,ie,n,nc,ii,io)
c********************************************************
c  Reads n data values from nc columns using * format
c   x - input data matrix
c========================================================
      real*8 x(n,nc)
      character mk
      intent(in) is,ie,n,nc,ii,io
      intent(out) x
c********************************************
      rewind(ii)
c  Skip header and parameter records
      read(ii,*)
      read(ii,*)
      if(is==1) then
       do i=1,n
        read(ii,*,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,*,end=1,err=2,iostat=ko)
       enddo
       do i=is,ie
        read(ii,*,end=1,err=2,iostat=ko)(x(i-is+1,k),k=1,nc)
       enddo
      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(9,*) '**** eof in read from data file '
      write(9,*) 'i = ',i
c  Error in data read
  2   write(io,*) '** Error in read of data file ',ko
      write(6,*) '** Error in read of data file ',ko
      write(9,*) '** Error in read of data file ',ko
      stop
      end
c
c===========================================
c
      subroutine wrh(su,n,p,pcu,pcl,pv,th,bpf,fl,fu,iunit,io)
c********************************************************
c  Writes run information for main output in output files
c========================================================
      integer p
      character*10 su
      character bpf
      intent(in) su,n,p,pcu,pcl,pv,th,bpf,fl,fu,iunit,io
c**********************
      if(iunit==1) then
c  Lower and upper frequencies for band
       if(bpf=='f') then
        write(io,'(/7x,''Passband ('',f9.3,2x,f10.3,'' ) kHz'')')fl,fu
       endif
      endif
c------
      if(iunit==2) then
       if(bpf=='f') then
        write(io,'(/7x,''Passband ('',f7.2,2x,f9.2,'' ) Hz'')')fl,fu
       endif
      endif
c------
      if(iunit==0) then
       if(bpf=='f') then
        if(fl<=1.e-17) then
         gl=float(n)
          else
           gl=fl
        endif
        if(fu<=1.e-17) then
         gu=2.
          else
           gu=fu
        endif
        write(io,'(/7x,''Passband ('',f9.2,2x,f9.2,'' )'',1x,a7)')gl,gu,
     *su
       endif
      endif
      if(pcu>0.0) then
       write(io,'(/4x,''The Right Tail of the Data is Trimmed'')')
       write(io,'(/4x,''Quantile to Trim Large Outliers = '',f4.1,''%'')
     *')pcu
      endif
      if(pcl>0.0) then
       write(io,'(/4x,''The Left Tail of the Data is Trimmed'')')
       write(io,'(/4x,''Quantile to Trim Small Outliers = '',f4.1,''%'')
     *')pcl
      endif
      if(p>0) then
       write(io,'(/4x,''Probability tail value threshold for the AR fit
     *iterations = '',f7.4/)')pv
      endif
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 = '',g15.3,7x,''Min value = '',g15.3)')max
     *,min
      write(io,'(60(''+''))')
c-------
      return
      end
c
c===========================================
c
      subroutine olsrun(y,x,r,n,p,ier,io)
c********************************************
c  Ordinary least squares fit of y=cx+e. Version 8-6-2009
c  x - array of nxp data matrix in column major order, y - dependent variable
c  c - estimates of c coefficients, am - estimate of intercept
c  r(1),...,r(n) - residuals, a - covariance matrix of estimates
c  r2 - Adjusted R square, se - stnd. deviation of residuals
c
c   ier = 1: A is not positive definite, ier = 2: p < 1,
c   ier = 3: p > n/2, ier = 4: work array can not be allocated
c============================================
      allocatable::a,b,c,d,m,v,w
      integer p
      real*8 x(n*p),y(n),r(n),m(:),v(:),w(:),c(:),a(:,:),b(:),d(:)
      intent(in) x,y,n,p,io
      intent(out) r,ier
c********************************************
c  Trap p
      if(p<1) then
       write(io,*)'p in OLS(p) ',p,' < 1'
       ier=2
       return
        elseif(p>n/2) then
         write(io,*)'p in OLS(p)',p,' > n/2 ',n/2,' n = ',n
         ier=3
         return
      endif
c-------
      na=p*p-1
      np=p-1
c  Allocate space
      allocate(a(0:np,0:np),b(0:np),d(0:np),m(0:np),c(0:np),v(0:na),w(0:
     *np),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays in olsrun'
       write(6,*)'Unable to allocate work space for arrays in olsrun'
       ier=4
       return
      endif
      call sls(y,x,c,a,r,v,b,w,d,m,n,p,r2,am,se,ier,io)
       write(io,'(/7x,''Data is Detrended using a Least Squares Fit of a
     * Centered Line''/)')
c  Write R square
      write(io,'(7x,''Adjusted R Square = '',f7.3,2x,''Standard Error of
     * Detrended Data = '',f10.4/)')r2,se
      write(io,'(7x,''Estimate of Intercept = '',f9.3/)')am
c  Write out estimates
      do k=0,np
       seb=sqrt(sngl(a(k,k)))
       tst=sngl(c(k))/seb
       write(io,'(2x,''b^('',i2,'') = '',f9.3,2x,''Sigma b^ = '',f12.4,2
     *x,''t = '',f10.2/)')k+1,c(k),seb,tst
      enddo
c-------
      deallocate(a,b,c,d,m,v,w)
      return
      end
c
c============================================
c
      subroutine sls(y,x,c,a,r,v,b,w,d,m,n,p,r2,am,se,ier,io)
c********************************************
c  Simple least squares (OLS) fit of y=bx+e
c============================================
      integer p,t
      real*8 y(n),x(n*p),m(p),r(n),v(p*p),w(p),c(p),a(p*p),b(p),d(p),dy,
     *sy,ssy,sse,sm
      intent(in) x,y,n,p,io
      intent(out) c,a,r,r2,am,se,ier
c********************************************
      call dcva(y,x,a,v,b,w,m,sy,ssy,n,p,io)
c  Cholesky decomposition of a
      call dcholdc(a,p,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       ier=1
       return
      endif
      call dcholsl(a,p,d,b,c)
c  Inverse of A
      call dinverse(a,p,d)
c-------
c  Residuals
      sse=0.d0
      do t=1,n
       sm=0.d0
       do k=1,p
        kt=(k-1)*n+t
        sm=sm+c(k)*dble(x(kt)-m(k))
       enddo
       dy=dble(y(t))-sy-sm
       sse=sse+dy*dy
       r(t)=sngl(dy)
      enddo
c  Standard error of residuals
      se=sqrt(sngl(sse)/float(n-p))
c  Intercept
      cn=0.
      do k=1,p
       cn=cn+sngl(c(k))*m(k)
      enddo
      am=sngl(sy)-cn
c  R squared
      r2=1.-sngl(sse/ssy)
      r2=r2-float(p)/float(n-p-1)*(1-r2)
c  Covariance matrix of estimates in a
      do j=1,p
       do k=j,p
        a((k-1)*p+j)=se*a((k-1)*p+j)
        a((j-1)*p+k)=a((k-1)*p+j)
       enddo
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine returns(y,r,n,io)
c********************************************
c  Rates of returns => r
c============================================
      allocatable::x2
      real*8 y(n),r(n),x2(:),yy
      integer t,n1
      intent(in) y,n,io
      intent(out) r
c*******************************************
c  Allocate space
      allocate(x2(0:n-1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays in returns'
       write(6,*)'Unable to allocate work space for arrays in returns'
       stop
      endif
c-------
c  Copy source data to scratch array x2
      do t=1,n
       x2(t-1)=0.d0
      enddo
      x2(0)=y(1)
      x2(n-1)=y(n)
      do t=1,n-2
       x2(t)=y(t+1)
      enddo
c  Calculate returns
      r(1)=0.d0
      do t=1,n-1
       yy=x2(t-1)
       if(yy<=0.d0) then
        r(t+1)=0.d0
         elseif((dabs(x2(t)/yy))<1.d-17) then
          write(io,'(2x,''return in Returns < 1.e-17 for t = '',i6)')t
          write(6,'(2x,''return in Returns < 1.e-17 for t = '',i6)')t
          write(9,'(2x,''return in Returns < 1.e-17 for t = '',i6)')t
          stop
           else
            r(t+1)=100.d0*dlog(x2(t)/yy)
       endif
      enddo
      deallocate(x2)
c-------
      return
      end
c
c===========================================
c
      subroutine difference(y,r,n,io)
c********************************************
c  Differences => r
c============================================
      implicit none
      real*8 y(n),r(n)
      integer n,t1,t,io
      intent(in) y,n
      intent(out) r
c********************************************
      do t=2,n
       t1=t-1
       r(t1)=y(t)-y(t1)
      enddo
      r(n)=0.d0
c-------
      return
      end
c
c===========================================
c
      subroutine hardclip(y,r,n,mean,io)
c********************************************
c  Hard clipped data => r
c============================================
      implicit none
      real*8 y(n),r(n),mean
      integer n,t,io
      intent(in) y,n,mean
      intent(out) r
c********************************************
      do t=1,n
       if(y(t)>=mean) then
        r(t)=1.d0
         else
          r(t)=-1.d0
       endif
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine dataclp(x,cth,n,io)
c********************************************
c     If x(t)<=cth then x(t)=x(t-1). If there is a string of x(t)'s<=cth
c     then x(t) is interpolated between the ends of the the string.
c     Version 8-7-2009
c============================================
      implicit none
      allocatable::x1,x2
      real*8 x(n),x1(:),x2(:),cth
      integer n,t,n1,ibad,io
      intent(in) n,cth,io
      intent(in out) x
c*******************************************
      n1=n-1
c     Allocate space
      allocate(x1(0:n1),x2(0:n1),stat=ibad)
      if(ibad/=0) then
       write(io,*)'Unable to allocate work space for arrays in dataaug'
       write(6,*)'Unable to allocate work space for arrays in dataaug'
       stop
      endif
c-------
c  Copy source data to x1
      do t=1,n
       x1(t-1)=x(t)
       x2(t-1)=cth
      enddo
c-------
      do t=2,n
       if(x(t)<=cth) then
        x(t)=x(t-1)  ! replace with previous value
       endif
      enddo
c-------
      x2(0)=x(1)
      x2(n-1)=x(n)
      do t=1,n-2
       x2(t)=x(t+1)
       if(x1(t)<=cth) then
        x2(t)=(x(t)+x(t+2))/2.d0  ! interpolate
       endif
      enddo
c  Clipped data => x
      do t=1,n
       x(t)=x2(t-1)
      enddo
c-------
      deallocate(x1,x2)
      return
      end
c
c===========================================
c
      subroutine dport(x,n,e,c,io)
c********************************************
c  c - portmentau test statistic p-value
c  x - data array, n - sample size, nlg=n**e, sum2 - sum of squared cors
c============================================================
      real*8 x(n),cv,sum2,rn
      real c,e
      integer t
      intent(in) x,n,e,io
      intent(out) c
c************************
      rn=dfloat(n-1)
      nlg=nint(float(n)**e)
      rng=float(nlg)
      c=1.
      sum2=0.d0
c  Start loop on k
      do k=1,nlg
       rn=dfloat(n-k)
c  Sum x(t)*x(t+k)
        cv=0.d0
        do t=1,n-k
         cv=cv+x(t)*x(t+k)
        enddo
       sum2=sum2+cv/rn*cv
      enddo
c  Chi**2 for squared cors
c  Correlation statistic
      call cdchi(sngl(sum2),rng,c,io)
      c=1.0-c
c--------
      return
      end