      program birrp
c
c  Revision 1 Oct 2001 added robust ar filter computation as default
c  Revision 2 Nov 2001 added basic and advanced input options, input of all components from one
c                      ascii file,use of dimensional period instead of Nyquists for input and
c                      output, output of apparent resistivity and phase
c  Revision 3 Mar 2002 numerous minor changes
c                      bug fixes in output section, provision of both period and frequency
c                      in output, allowance of sample rate as well as sample interval,
c                      change to scale estimate for Thomson weights
c  Revision 3.1        fixed minor errors in rarfilt
c  Revision 3.2        changed convergence flag to integer
c  Revision 3.3        added immediate termination if all weights go to zero
c  Revision 4 Dec 2002 major improvements to robust ar filter routine, including exclusive
c                      use of double precision to get acvs and ar filter
c                      ar filter is checked for minimum phase behavior
c                      ar filter and poles output in diagnostic file
c                      added option to input prewhitened data and prewhitening filter
c  Revision 5 Mar 2004 Major code revision
c                      changed convention for frequency input to FILTER to dimensional units, added
c	               automatic interpolation feature from filter data input from file when nfil<0
c                      added J-format output, added option for automatic computation of data bounds
c                      given start and end times of each time series, integrated selection of
c                      delete one files with output level, added sort feature for responses in
c                      case frequencies are not increasing, added diagnostic output for first
c                      stage, made minor corrections to qq plot output, removed l1 initial solution
c                      option, substantially reduced memory usage, added full output of transfer function,
c                      jackknife errors, and coherences for stage 1
c  Revision 5.1        Fixed coherence thresholding so frequency selection works
c                      forced fft size in call to rarfilt to be even to 
c                      make levinson recursion stable
c  Revision 5.2        Changed to Bonferroni confidence intervals for
c                      apparent resistivity and phase
c  Revision 5.2.1      Corrected fatal error in nonrobust processing using 2 stage approach
c  Revision 5.2.2      Corrected indexing error in ccoher1 and ccoher2
c  Revision 5.2.3      Changed zeromean of data segments to zero median
c                      Changed to direct median computation in mad
c
c  Computes transfer function and multiple coherence between up to
c  three output and three input data series using bounded influence,
c  two stage least squares method.
c
c  Computation is by section averaging using variable section lengths.
c  Robust weights are computed using a Rayleigh distribution model on the
c  regression residuals and applied equally to the real and imaginary parts.
c  Leverage weights are computed by comparing the hat matrix diagonal elements 
c  to their expected value. Optional coherence thresholding can be applied with
c  a combination of a selected coherence value and/or a fractional trim value.
c  Confidence limits are obtained using the unbalanced jackknife.
c
c  parameter list in 'parameters.h' appears below
c
c  inter = logical flag for interactive (.true.) or runfile (.false.) version
c  npcsm = maximum number of data pieces
c  nptsm = maximum number of data in all pieces
c  nptssm = maximum fft section length
c  noutm = maximum number of output time series
c  ninpm = maximum number of input time series
c  nrefm = maximum number of reference time series, must be .ge. ninpm
c  nrsitem = maximum number of reference sites
c  nsectm = maximum number of data sections, should be set to
c           nptsm*nsctinc**(nsctmax-1)/(nptssm*offs), where
c           offs=0.57-0.07*tbw and tbw is the time bandwidth
c
      include 'parameters.h'
      parameter (nserm=noutm+ninpm+nrefm)
c
c  additional parameters which should not need to be changed frequently
c
c  nblk = block size for coherence thresholding, set to 10
c            if changed, it must also be adjusted in bcoher2_e and bcoher2_b
      parameter (nblk=10)
c  nfsm = maximum number of frequencies to be used for a given section size
      parameter (nfsm=2)
c  nfm = maximum number of output frequencies
      parameter (nfm=100)
c  nfilm = maximum number of filter parameters
      parameter (nfilm=100)
c         if changed, it must also be adjusted in dataft
c  narm = maximum ar filter length
c         if changed, it must also be adjusted in dataft, rarfilt1, and binput
      parameter (narm=100)
c
c  mfftm = maximum value of mfft
      parameter (mfftm=10)
c
c  nrecl = maximum length of a binary record, set to 30000 in subroutine input
c
      character filnam*80,fildat(nserm,npcsm)*80,frm*1,ofil*70,
     $          postfix*30,cn*3,cpar(nserm)*80,jlabel1(2,2)*3,
     $          jlabel2(3,2)*8,dstim(nserm,npcsm)*19,
     $          wstim(nserm,npcsm)*19,wetim(nserm,npcsm)*19,
     $          cvar(3)*20,cvar1(3)*20
      integer tdiff,dssec,wssec,wesec,wsmax,wemin
      dimension dssec(nserm,npcsm),wssec(nserm,npcsm),
     $          wesec(nserm,npcsm),wsmax(npcsm),wemin(npcsm)
      double complex b,aa,bb,ob,root
      dimension b(nrefm),ob(nrefm),aa(nrefm,nrefm),bb(nrefm),
     $          root(narm+1)
      complex fft,ffts,tbn,tbr,hn,hr,hd1,gn1,gn2,gr1,gr2,resq,oresq
      dimension fft(nptssm),ffts(nsectm,nserm,nfsm),
     $          tbn(ninpm,nrefm,nfm),hn(noutm,ninpm,nfm),
     $          tbr(ninpm,nrefm,nfm),hr(noutm,ninpm,nfm),
     $          gn1(ninpm,nfm),gn2(noutm,nfm),
     $          gr1(ninpm,nfm),gr2(noutm,nfm),
     $          gn1lo(ninpm,nfm),gn1hi(ninpm,nfm),
     $          gr1lo(ninpm,nfm),gr1hi(ninpm,nfm),
     $          gn2lo(noutm,nfm),gn2hi(noutm,nfm),
     $          gr2lo(noutm,nfm),gr2hi(noutm,nfm),
     $          dfn1(ninpm,nfm),dfn2(noutm,nfm),
     $          dfr1(ninpm,nfm),dfr2(noutm,nfm),
     $          hd1(nrefm+1,nsectm),gd1(nsectm)
      double precision sum,var,wtsum,dpss,dwork,dofscal,evalu,coef
      dimension dpss(nptssm),dwork(6*nptssm+13),coef(narm+1)
      dimension frej(nfm),rdata(nptsm,nserm),
     $          work(14*nptsm+7*nptsm*mfftm/nptssm+18*nptssm+32),
     $          fpar(nfilm+1,nserm),arf(narm,nserm,npcsm),
     $          wsave(4*nptssm+15),hat(nsectm),
     $          res(nsectm),pwro(nsectm),pwri(nsectm),
     $          xindex(nsectm),w(nsectm),
     $          what(nsectm),wt(nsectm),wt1(nsectm)
      dimension varn1(ninpm,nrefm,nfm),varn2(noutm,ninpm,nfm),
     $          varr1(ninpm,nrefm,nfm),varr2(noutm,ninpm,nfm)
      dimension temp(nsectm),freq(nfm),arft(narm),
     $          theta1(nrsitem+2),theta2(nrsitem+2),
     $          phi(nrsitem+2),c2b(nsectm/nblk,ninpm),
     $          c2e(nsectm/nblk,noutm),
     $          fracthresh1(ninpm,nfm),fracdat1(ninpm,nfm),
     $          fracthresh2(noutm,nfm),fracdat2(noutm,nfm)
      integer convflag1,convflag2
      dimension nstart(npcsm+1),npts(npcsm),nfil(nserm),
     $          ncomp(nserm),iindex(noutm+ninpm,nserm),
     $          nskip(nserm,npcsm),nread(npcsm),nblock(nserm),
     $          nfr(nfm),convflag1(ninpm,nfm),convflag2(noutm,nfm)
      logical compare,c2flag,c2flag1,cflag1,cflag2,
     $        lwt,lwt1,llo,lhi,llo1,lhi1,ilev
      dimension c2flag(nsectm/nblk),c2flag1(nsectm/nblk),
     $          lwt(nsectm),lwt1(nsectm),llo(nsectm),lhi(nsectm),
     $          llo1(nsectm),lhi1(nsectm)
      double precision pi,rad
      parameter (pi=3.141592653589793d0,rad=pi/180.d0)
      equivalence (ffts,dwork,work)
      data nrr/0/,nrr1/0/,nr2/1/,nr3/0/,c2threshe/0./,c2threshe1/0./,
     $     c2threshb/0./,nz/0/,nblock/nserm*0/,mfft/2/,nprej/0/,
     $     cn/'   '/
      data f2lo/0./,f2hi/0.5/,eps/0.01/,eps1/0.1/
      data jlabel1(1,1)/'RXX'/,jlabel1(1,2)/'RXY'/,
     $     jlabel1(2,1)/'RYX'/,jlabel1(2,2)/'RYY'/
      data jlabel2(1,1)/'ZXX S.I.'/,jlabel2(1,2)/'ZXY S.I.'/,
     $     jlabel2(2,1)/'ZYX S.I.'/,jlabel2(2,2)/'ZYY S.I.'/,
     $     jlabel2(3,1)/'TZX     '/,jlabel2(3,2)/'TZY     '/
      data ilev/.false./
c
c  input run parameters
c
      if(inter)write(6,5)
      read(5,*)nlev
      if(nlev.eq.1)ilev=.true.
      if(inter)write(6,6)
      read(5,*)nout
      if(inter)write(6,7)
      read(5,*)ninp
      if(ilev)then
        if(inter)write(6,8)
        read(5,*)nref
        if(nref.gt.3)then
          nrr=1
          if(inter)write(6,9)
          read(5,*)nr3,nr2
        elseif(nref.eq.3)then
          if(ninp.lt.3)nrr=1
          nr3=1
          nr2=0
        endif
      else
        nref=ninp
      endif
      nser=nout+ninp+nref
      if(ilev)then
        if(inter)write(6,10)
        read(5,*)nrr1
        if((nrr.eq.1).and.(nrr1.eq.0))then
          if(inter)write(6,*)' 2 stage processing required with this',
     $                       ' input/output combination--set to 1'
        nrr1=1
        endif
        nrr=nrr1 
      endif
    1 if(inter)write(6,11)
      read(5,*)tbw
      if(inter)then
        if((tbw.lt.1.).or.(tbw.gt.4.))then
          write(6,*)' time-bandwidth out of bounds-should be 1 to 4'
          goto 1
        endif
      endif
      offs=0.57-0.07*tbw
      if(inter)write(6,12)
      read(5,*)deltat
      if(deltat.eq.0.)deltat=1.
      if(deltat.lt.0.)deltat=1./abs(deltat)
    2 continue
      if(ilev)then
        if(inter)write(6,13)
        read(5,*)nfft,nsctinc,nsctmax
      else
        if(inter)write(6,14)
        read(5,*)nfft,nsctmax
        nsctinc=2
      endif
      if(mod(nfft,2).ne.0)nfft=nfft+1
      if(ilev)then
        if(inter)write(6,15)
        read(5,*)nf1,nfinc,nfsect
      else
        nf1=nint(tbw+2.)
        nfinc=nint(tbw)
        nfsect=2
      endif
      if(inter)then
        if(nptsm*float(nsctinc**(nsctmax-1))/(nfft*offs).gt.nsectm)then
          write(6,*)' number of sections exceeds dimension of nsectm'
          goto 2
        endif
        if(nfft.gt.nptssm)then
          write(6,*)' fft length exceeds dimension of nptssm'
          goto 2
        endif
        if(nsctmax*nfsect.gt.nfm)then
          write(6,*)' number of periods exceeds dimension of',nfm
          goto 2
        endif
        if(nfsect.gt.nfsm)then
          write(6,*)' periods per section exceeds dimension of',nfsm
          goto 2
        endif
        write(6,*)' table of output periods by section length'
        nfft1=nfft
        do 3 i=1,nsctmax
          write(6,*)(deltat*nfft1/(nf1+j*nfinc),j=0,nfsect-1)
   3      nfft1=nfft1/nsctinc
        write(6,16)
        read(5,'(a)')frm
        if((frm.eq.'n').or.(frm.eq.'N'))goto 2
      endif
      if(ilev)then
        write(6,17)
        read(5,*)mfft
      endif
      if(nfft/mfft.lt.nptssm/mfftm)then
        mfft=nfft*mfftm/nptssm
        if(mfft.eq.0)then
          write(6,*)' mfftm is too small for input value of nfft'
          stop
        endif
        write(6,*)' mfft adjusted to',mfft
      endif 
      if(ilev)then
        if(inter)write(6,18)
        read(5,*)uin,ainlin,ainuin
      else
        ainlin=-999
        write(6,19)
        read(5,*)uin,ainuin
      endif
      if(nrr.eq.1)then
        if(inter)write(6,20)
        read(5,*)c2threshb
        c2threshb=min(c2threshb,1.)
      endif
      if(inter)write(6,21)
      read(5,*)c2threshe
      c2threshe=min(c2threshe,1.)
      if(nout.eq.3)then
        if(inter)write(6,22)
        read(5,*)nz
        if(nz.eq.0)then
          if(inter)write(6,21)
          read(5,*)c2threshe1
        endif
      endif
      if(ilev)then
        if((c2threshb.ne.0.).or.(c2threshe.ne.0.).or.
     $     (c2threshe1.ne.0))then
           if(inter)write(6,23)
          read(5,*)perlo,perhi
          if(perlo.lt.0)perlo=1./abs(perlo)
          if(perhi.lt.0.)perhi=1./abs(perhi)
          f2lo=deltat/perlo
          f2hi=min(0.5,deltat/perhi)
        endif
      endif
      c2crit=1.-0.05**(2./(2.*(nblk-ninp)))
      if((inter).and.(c2threshb.ne.0.).and.(c2threshb.lt.c2crit))then
        write(6,*)' zero coherence at 95% significance of',c2crit
        write(6,*)' is larger than specified value of ',c2threshb
        write(6,*)' for 1st stage'
      endif
      if((inter).and.(c2threshe.ne.0.).and.(c2threshe.lt.c2crit))then
        write(6,*)' zero coherence at 95% significance of',c2crit
        write(6,*)' is larger than specified value of ',c2threshe
        write(6,*)' for 2nd stage'
      endif
      if((inter).and.(nz.eq.0).and.(c2threshe1.ne.0.).and.
     $   (c2threshe1.lt.c2crit))then
        write(6,*)' zero coherence at 95% significance of',c2crit
        write(6,*)' is larger than specified value of ',c2threshe1
        write(6,*)' for third output variable'
      endif
      if(inter)write(6,24)
      read(5,'(a)')ofil
      call strlen(ofil,l1,l2)
c
c  open diagnostic file
c
      open (unit=19,file=ofil(l1:l2)//'.diag',
     $      status='unknown')
      if(.not.inter)then
        if((tbw.lt.1.).or.(tbw.gt.4.))then
          write(19,*)' ERROR-time-bandwidth out of bounds-',
     $                 ' should be 1 to 4'
          stop
        endif
        if(nptsm*float(nsctinc**(nsctmax-1))/(nfft*offs).gt.nsectm)then
          write(19,*)' ERROR-number of sections exceeds',
     $                 ' dimension of',nsectm
          stop
        endif
        if(nfft.gt.nptssm)then
          write(19,*)' ERROR-fft length exceeds dimension',
     $                 ' of',nptssm
          stop
        endif
        if(nsctmax*nfsect.gt.nfm)then
          write(19,*)' ERROR-number of frequencies exceeds',
     $                 ' dimension of',nfm
          stop
        endif
        if(nfsect.gt.nfsm)then
          write(19,*)' ERROR-frequencies per section exceeds',
     $                 ' dimension of',nfsm
          stop
        endif
      endif
      if(inter)write(6,25)
      read(5,*)nlev
      if(ilev)then
        if(inter)write(6,26)
        read(5,*)nprej
        if(nprej.ne.0)then
          if(inter)write(6,27)
          do 4 i=1,nprej
            read(5,*)prej
            if(prej.lt.0.)prej=1./abs(prej)
    4       frej(i)=deltat/prej
        endif
      endif
    5 format(' select the input level (0=basic,1=advanced):',$)
    6 format(' input the number of output time series:',$)
    7 format(' input the number of input time series:',$)
    8 format(' input the number of reference time series:',$)
    9 format(' input the number of 3/2 component reference sites:',$)
   10 format(' select r-r (0) or 2-stage robust method (1):',$)
   11 format(' select tbw for prolate data window (1 to 4):',$)
   12 format(' input the data sample interval:',$)
   13 format(' input the initial section length, divisor,',
     $       ' and maximum number of sections:',$)
   14 format(' input the initial section length and',
     $       ' maximum number of sections:',$)
   15 format(' input the starting index, increment, and number of',
     $       ' frequencies per section:',$)
   16 format(' are these values acceptable (y or n)?',$)
   17 format(' input the robust ar filter section length divisor:',$)
   18 format(' input the robustness and lower/upper leverage',
     $       ' parameters (-999 to omit):',$)
   19 format(' input the robustness and leverage parameters',
     $       ' (0 and >0.99 are recommended):',$)
   20 format(' input the 1st stage coherence threshold:',$)
   21 format(' input the 2nd stage coherence threshold:',$)
   22 format(' input the 3rd variable threshold mode',
     $       ' (0 to threshold separately, 1 to threshold using',
     $       '  2nd stage value, and 2 to threshold simultaneously',
     $       '  with 2nd stage variables):',$)
   23 format(' input the longest and shortest periods:',$)
   24 format(' input the output filename root:',$)
   25 format(' input the output level (-3 to 3):',$)
   26 format(' input the number of rejection frequencies:',$)
   27 format(' input the frequencies')
c
c  input data
c
      if(inter)write(6,50)
      read(5,*)npcs
      if(inter)write(6,51)
      read(5,*)nar
      if(inter)write(6,52)
      read(5,*)imode
      if(inter)write(6,53)
      read(5,*)jmode
      if(imode.eq.0)then
        frm='f'
      elseif(imode.eq.1)then
        frm='u'
      endif
      if(uin.ne.-999.)then
c output header block for j-file
        filnam=ofil(l1:l2)//'.j'
	open(unit=18,file=filnam,status='unknown')
	call jfinit(ilev,nout,ninp,nref,nr3,nr2,nrr,
     $              tbw,deltat,nfft,nsctinc,nsctmax,nf1,
     $              nfinc,nfsect,mfft,uin,ainlin,ainuin,
     $              c2threshb,c2threshe,nz,c2threshe1,
     $              perlo,perhi,nprej,frej,npcs,nar,imode,
     $              jmode)
      endif
c
c input data parameters
c
      if(imode.lt.2)then
        do 36 nn=1,npcs
	  if(jmode.eq.0)then
	    if(inter)write(6,54)
	    read(5,*)nread(nn)
	  endif
          do 36 n=1,nser
            if(nn.eq.1)then
              if(inter)write(6,55)
              read(5,*)nfil(n)
              if(nfil(n).gt.nfilm)then
                if(inter)then
                  write(6,*)' number of filter coefficients exceeds',
     $                      ' dimension'
                else
                  write(19,*)' ERROR-number of filter coefficients',
     $                         ' exceeds dimension'
                endif 
                stop
              endif
              if(nfil(n).gt.0)then
                if(inter)write(6,56)
                read(5,*)(fpar(i,n),i=1,nfil(n))
	        fpar(nfil(n)+1,n)=deltat
              elseif(nfil(n).lt.0)then
	        if(inter)write(6,57)
	        read(5,'(a)')cpar(n)
	      endif
              if(imode.eq.1)then
                if(inter)write(6,58)
                read(5,*)nblock(n)
              endif
	    endif
            if(nar.lt.0)then
              if(inter)write(6,59)
              read(5,'(a)')filnam
              open(unit=10,file=filnam,status='old')
              do 33 i=1,-nar
   33           read(10,*,end=34)arf(i,n,nn)
              close(unit=10)
              goto 35
   34         if(inter)then
                write(6,*)' end of file for ar filter at series',n,
     $                    ' piece',nn
              else
                write(19,*)' ERROR-end of file for ar filter at',
     $                       ' series',n,' piece',nn
              endif
              stop
   35         continue
            endif
            if(inter)write(6,60)
            read(5,'(a)')fildat(n,nn)
            if(jmode.eq.0)then
              if(inter)write(6,61)
              read(5,*)nskip(n,nn)
            else
              if(inter)write(6,62)
              read(5,'(a19)')dstim(n,nn)
              read(5,'(a19)')wstim(n,nn)
              read(5,'(a19)')wetim(n,nn)
            endif
   36       continue
      else
        do 37 n=1,nser
          if(inter)write(6,55)
          read(5,*)nfil(n)
          if(nfil(n).gt.nfilm)then
            if(inter)then
              write(6,*)' number of filter coefficients exceeds',
     $                  ' maximum dimension'
            else
              write(19,*)' ERROR-number of filter coefficients',
     $                     ' exceeds maximum dimension'
            endif 
            stop
          endif
          if(nfil(n).gt.0)then
            if(inter)write(6,56)
            read(5,*)(fpar(i,n),i=1,nfil(n))
	    fpar(nfil(n)+1,n)=deltat
          elseif(nfil(n).lt.0)then
	    if(inter)write(6,57)
	    read(5,'(a)')cpar(n)
	  endif
   37     continue
        if(inter)write(6,63)
        read(5,*)ncomp(1),(iindex(i,1),i=1,nout+ninp)
        do 38 nn=1,npcs
          if(inter)write(6,60)
          read(5,'(a)')fildat(1,nn)
          if(jmode.eq.0)then
            if(inter)write(6,61)
            read(5,*)nskip(1,nn)
            if(inter)write(6,54)
            read(5,*)nread(nn)
          else
            if(inter)write(6,62)
            read(5,'(a19)')dstim(1,nn)
            read(5,'(a19)')wstim(1,nn)
            read(5,'(a19)')wetim(1,nn)
          endif
   38     continue
        do 41 n=1,nr3
          if(inter)write(6,63)
          read(5,*)ncomp(n+1),(iindex(i,n+1),i=1,3)
          do 41 nn=1,npcs
            if(inter)write(6,60)
            read(5,'(a)')fildat(n+1,nn)
            if(jmode.eq.0)then
              if(inter)write(6,61)
              read(5,*)nskip(n+1,nn)
            else
              if(inter)write(6,62)
              read(5,'(a19)')dstim(n+1,nn)
              read(5,'(a19)')wstim(n+1,nn)
              read(5,'(a19)')wetim(n+1,nn)
            endif
   41       continue
        do 43 n=1,nr2
          if(inter)write(6,63)
          read(5,*)ncomp(n+nr3+1),(iindex(i,n+nr3+1),i=1,2)
          do 43 nn=1,npcs
            if(inter)write(6,60)
            read(5,'(a)')fildat(n+nr3+1,nn)
            if(jmode.eq.0)then
              if(inter)write(6,61)
              read(5,*)nskip(n+nr3+1,nn)
            else
              if(inter)write(6,62)
              read(5,'(a19)')dstim(n+nr3+1,nn)
              read(5,'(a19)')wstim(n+nr3+1,nn)
              read(5,'(a19)')wetim(n+nr3+1,nn)
            endif
   43       continue
      endif
c
   50 format(' input the number of data pieces:',$)
   51 format(' input the length of the ar filter (0 for none,'
     $       ' <0 for filename):',$)
   52 format(' select the file mode (0=separate ascii files,',
     $       ' 1=separate binary files',
     $       ' 2=all component ascii files',
     $       ' 3=ts ascii file:',$)
   53 format(' select the input mode (0=by points, 1=by date/time):',$)
   54 format(' input the number of points to read:',$)
   55 format(' input the number of filter parameters',
     $       ' (<0 for filename):',$)
   56 format(' input the filter parameters')
   57 format(' input the filter filename:',$)
   58 format(' input the block size:',$)
   59 format(' input the ar filter filename:',$)
   60 format(' input the data filename:',$)
   61 format(' input the number of points to skip:',$)
   62 format(' input the data start and window start/end times'
     $       ' (yyyy-mm-dd hh:mm:ss):',$)
   63 format(' input the number of components and the indices',
     $       ' by component:',$)
c
c align time series if date/times were input
c
      if(jmode.eq.1)then
c
c  reference time differences to 1980-01-01 00:00:00
c
        if(imode.lt.2)then
          m=nser
        else
          m=1+nr3+nr2
        endif
        do 65 nn=1,npcs
         wsmax(nn)=0
         wemin(nn)=86400*365*50
         do 65 n=1,m
            dssec(n,nn)=tdiff('1980-01-01 00:00:00',dstim(n,nn))
            wssec(n,nn)=tdiff('1980-01-01 00:00:00',wstim(n,nn))
            wesec(n,nn)=tdiff('1980-01-01 00:00:00',wetim(n,nn))
            wsmax(nn)=max(wsmax(nn),wssec(n,nn))
   65       wemin(nn)=min(wemin(nn),wesec(n,nn))
        do 66 nn=1,npcs
          nread(nn)=dble(wemin(nn)-wsmax(nn))/dble(deltat)
          if(nread(nn).le.0)then
            if(inter)then
              write(6,*)' nread not positive for piece',nn
            else
              write(19,*)' nread not positive for piece',nn
            endif
            stop
          endif
          do 66 n=1,m
   66       nskip(n,nn)=dble(wsmax(nn)-dssec(n,nn))/dble(deltat)
      endif
c
c  input data
c
      nstart(1)=1
      nptstot=0
      if(imode.lt.2)then
        do 76 n=1,nser
          if(uin.ne.-999.)then
            write(18,996)'#nfil=',nfil(n)
          endif
          do 76 nn=1,npcs
            if(uin.ne.-999.)then
               write(18,'(2a)')'#filnam=',fildat(n,nn)
              write(18,996)'#nskip=',nskip(n,nn),' nread=',nread(nn)
            endif
            call input(fildat(n,nn),frm,nblock(n),nskip(n,nn),
     $                 nread(nn),rdata(nstart(nn),n))
            if(nread(nn).eq.0)then
              if(inter)then
                write(6,*)' end of file at series',n,' piece',nn
              else
                write(19,*)' ERROR-end of file at series',n,
     $                       ' piece',nn
              endif
              stop
            endif
            if(n.eq.1)then
              nptstot=nptstot+nread(nn)
              npts(nn)=nread(nn)
              nstart(nn+1)=nstart(nn)+nread(nn)
              if(nptstot.gt.nptsm)then
                if(inter)then
                  write(6,*)' number of data exceeds dimension of',
     $                      nptsm
                else
                  write(19,*)' ERROR-number of data exceeds',
     $                         ' dimension of',nptsm
                endif
                stop
              endif
            elseif(nread(nn).ne.npts(nn))then
              if(inter)then
                write(6,*)' number of data incorrect at series',n,
     $                    ' piece',nn,' nread=',nread(nn)
              else
                write(19,*)' ERROR-number of data incorrect at',
     $                       ' series', n,' piece',nn,
     $                       ' nread=',nread(nn)
              endif
              stop
            endif
   76       continue
      else
        do 77 n=1,nser
          if(uin.ne.-999.)then
            write(18,996)'#nfil=',nfil(n)
          endif
   77     continue
        if(uin.ne.-999.)then
          write(18,996)'#ncomp=',ncomp(1)
          write(18,997)'#indices=',(iindex(i,1),i=1,nout+ninp)
        endif
        do 78 nn=1,npcs
          if(uin.ne.-999.)then
            write(18,'(2a)')'#filnam=',fildat(1,nn)
            write(18,996)'#nskip=',nskip(1,nn),' nread=',nread(nn)
          endif
          call binput(fildat(1,nn),imode-1,nskip(1,nn),nread(nn),
     $                nar,ncomp,nout+ninp,iindex(1,1),work,
     $                rdata(nstart(nn),1),arf(1,1,nn),cvar)
          if(nread(nn).eq.0)then
            if(inter)then
              write(6,*)' end of file at piece',nn
            else
              write(19,*)' ERROR-end of file at piece',nn
            endif
            stop
          endif
          nptstot=nptstot+nread(nn)
          if(nptstot.gt.nptsm)then
            if(inter)then
              write(6,*)' number of data exceeds dimension of',nptsm,
     $                  ' at piece',nn
            else
              write(19,*)' ERROR-number of data exceeds',
     $                     ' dimension of',nptsm,' at piece',nn
            endif
            stop
          endif
          npts(nn)=nread(nn)
   78     nstart(nn+1)=nstart(nn)+nread(nn)
        mser=nout+ninp+1
        do 82 n=1,nr3
          if(uin.ne.-999.)then
            write(18,996)'#ncomp=',ncomp(n+1)
            write(18,997)'#indices=',(iindex(i,n+1),i=1,3)
          endif
          do 81 nn=1,npcs
            if(uin.ne.-999.)then
              write(18,'(2a)')'#filnam=',fildat(n+1,nn)
              write(18,996)'#nskip=',nskip(n+1,nn),
     $                     ' nread=',nread(nn)
            endif
            call binput(fildat(n+1,nn),imode-1,nskip(n+1,nn),
     $                  nread(nn),nar,ncomp(n+1),3,iindex(1,n+1),work,
     $                  rdata(nstart(nn),mser),arf(1,mser,nn),cvar1)
             if(nread(nn).eq.0)then
              if(inter)then
                write(6,*)' end of file for reference',n,
     $                    ' at piece',nn
              else
                write(19,*)' ERROR-end of file for reference',n,
     $                       ' at piece',nn
              endif
              stop
            elseif(nread(nn).ne.npts(nn))then
              if(inter)then
                write(6,*)' number of data incorrect for reference',n,
     $                    ' at piece',nn
              else
                write(19,*)' ERROR-number of data incorrect for',
     $                       ' reference',n,' at piece',nn
              endif
              stop
            endif
   81       continue
   82     mser=mser+3
        do 84 n=1,nr2
          if(uin.ne.-999.)then
            write(18,996)'#ncomp=',ncomp(n+nr3+1)
            write(18,997)'#indices=',(iindex(i,n+nr3+1),i=1,2)
          endif
          do 83 nn=1,npcs
            if(uin.ne.-999.)then
              write(18,'(2a)')'#filnam=',fildat(n+nr3+1,nn)
              write(18,996)'#nskip=',nskip(n+nr3+1,nn),
     $                     ' nread=',nread(nn)
            endif
            call binput(fildat(n+nr3+1,nn),imode-1,
     $                  nskip(n+nr3+1,nn),nread(nn),nar,
     $                  ncomp(n+nr3+1),2,iindex(1,n+nr3+1),
     $                  work,rdata(nstart(nn),mser),
     $                  arf(1,mser,nn),cvar1)
            if(nread(nn).eq.0)then
              if(inter)then
                write(6,*)' end of file for reference',n+nr3+1,
     $                    ' at piece',nn
              else
                write(19,*)' ERROR-end of file for reference',
     $                       n+nr3+1,' at piece',nn
              endif
              stop
            elseif(nread(nn).ne.npts(nn))then
              if(inter)then
                write(6,*)' number of data incorrect for reference',
     $                    n+nr3+1,' at piece',nn
              else
                write(19,*)' ERROR-number of data incorrect for',
     $                       ' reference',n+nr3+1,' at piece',nn
              endif
              stop
            endif
   83       continue
   84     mser=mser+2
      endif
c
c  force data to have zero median for each piece unless data have been prewhitened
c
      if(nar.ge.0)then
        do 88 n=1,nser
          do 88 nn=1,npcs
            xmed=xmedian(rdata(nstart(nn),n),npts(nn))
            do 87 i=nstart(nn),nstart(nn)+npts(nn)-1
   87         rdata(i,n)=rdata(i,n)-xmed 
   88       continue
      endif
c
c  input rotation parameters
c
      if(inter)write(6,95)
      if(nout.ge.2)read(5,*)theta1(1),theta2(1),phi(1)
      if(uin.ne.-999.)write(18,998)'#theta1=',theta1(1),
     $                             ' theta2=',theta2(1),
     $                             ' phi=',phi(1)
      read(5,*)theta1(2),theta2(2),phi(2)
      if(uin.ne.-999.)write(18,998)'#theta1=',theta1(2),
     $                             ' theta2=',theta2(2),
     $                             ' phi=',phi(2)
      do 90 i=1,nr3+nr2
        read(5,*)theta1(2+i),theta2(2+i),phi(2+i)
   90   if(uin.ne.-999.)write(18,998)'#theta1=',theta1(2+i),
     $                               ' theta2=',theta2(2+i),
     $                               ' phi=',phi(2+i)
   95 format(' input the rotation angles (cw=+)')
c
c  open robust delete-one files
c
      if(nlev.lt.0)then
        do 162 nn=1,nout
          filnam=ofil(l1:l2)//'.'//char(48+nn)//'r.d1rf'
          open (unit=32+nn,file=filnam,form='unformatted',
     $          status='unknown')
  162   continue
      endif
c
c  prewhiten data series with robust ar filter if required
c
        if(nar.gt.0)then
c          nfftmin=nfft/(nsctinc**(nsctmax-1))
c  use average of logs of shortest and longest fft for ar filter ffts
c          nfft1=exp((log(float(nfft))+log(float(nfftmin)))/2.)
          nfft1=nfft/mfft
          if(mod(nfft1,2).ne.0)nfft1=nfft1+1
          write(19,*)' ar filter fft length=',nfft1
          do 164 jj=1,npcs
            n1=nstart(jj)
            do 163 j=1,nser
              n=npts(jj)
              call rarfilt(n,nfft1,nar,rdata(n1,j),work,arft,info)
              if(info.gt.0)then
                write(19,*)' WARNING-error in call to rarfilt',
     $                       ' series=',j,' piece=',jj,
     $                       ' info=',info
                if(inter)then
                  write(6,*)' Error in call to rarfilt',
     $                      ' series=',j,' piece=',jj,
     $                      ' info=',info
                endif
              endif
              call prewhiten(n,nar,rdata(n1,j),arft,rdata(n1,j))
              do 163 i=1,nar
  163           arf(i,j,jj)=arft(i)
  164       if(j.eq.nser)npts(jj)=n
        endif
c
c check ar filter for minimum phase behavior
c
        do 167 jj=1,npcs
          do 167 j=1,nser
            coef(1)=1.
            do 165 i=2,nar+1
  165         coef(i)=arf(i-1,j,jj)
            call dpqr79(nar,coef,root,ierr,work)
            write(19,*)' poles of ar z-transform and distance from',
     $                   ' unit circle for series=',j,' piece=',jj
            do 166 i=1,nar
              if(abs(root(i)).ge.1.)then
                if(inter)write(6,*)' Non-minimum phase ar filter for',
     $                             ' series ',j,' piece ',jj
                write(19,*)' WARNING-non minimum phase ar filter'
              endif
              x1=dreal(root(i))
              y1=dimag(root(i))
              if(x1.eq.0.)then
                dist=abs(1.-y1)
              else
                x0=sign(sqrt(x1*x1/(x1*x1+y1*y1)),x1)
                y0=y1*x0/x1
                dist=sqrt((x1-x0)**2+(y1-y0)**2)
              endif 
              write(19,*)x1,y1,dist
  166         continue
            write(19,*)' ar filter coefficients for series=',j,
     $                   ' piece=',jj
            write(19,*)(arf(i,j,jj),i=1,nar)
  167       continue
c
c  begin outer loop pivoting on section length
c
      nfl=1
      do 171 j=1,nfm
	do 170 n=1,nout
  170     convflag2(n,j)=0
        do 171 n=1,ninp
  171     convflag1(n,j)=0
      do 790 nsc=1,nsctmax
        write(19,*)'*******************************************'
        write(19,*)' section length=',nfft
        xnfft=nfft
c
c  initialize fft computation and compute slepian sequence data window
c
        call rtpss(nptssm,nfft,1,dble(tbw/xnfft),dpss,evalu,
     $             info,dwork)
        if(info.ne.0)then
          if(inter)write(6,*)' error in call to rtpss info=',info
          write(19,*)' WARNING-error in call to rtpss info=',info
        endif
        call cffti(nfft,wsave)
c
c  compute section ffts and correct for instrument response and prewhitening
c
        do 189 k=1,nfsm
          do 189 j=1,nserm
            do 189 i=1,nsectm
              ffts(i,j,k)=0.
  189         continue
        nsect=0
        noff=offs*xnfft+1
        do 190 jj=1,npcs
          if(npts(jj).le.nfft)goto 190
          nsct=float(npts(jj)-nfft-1)/float(noff)+1
          if(npts(jj)-((nsct-1)*noff+nfft+1).ge.nfft)nsct=nsct+1
          ko=nstart(jj)
          do 180 j=1,nser
            call dataft(nsct,j,nf1,nfinc,nfsect,nfft,noff,
     $                  nfil(j),abs(nar),deltat,dpss,wsave,
     $                  rdata(ko,j),fpar(1,j),cpar(j),arf(1,j,jj),
     $                  fft,ffts(nsect+1,1,1))
  180       continue
          nsect=nsect+nsct
  190     continue
        nsect=(nsect/nblk)*nblk
        if(nsect.gt.nsectm)then
          if(inter)then
            write(6,*)' number of sections ',nsect,
     $                ' is larger than maximum dimension of ',nsectm
          else
            write(19,*)' ERROR-number of sections ',nsect,
     $                   ' is larger than maximum dimension of ',nsectm
          endif
          stop
        elseif(nsect.eq.0)then
          goto 790
        endif
        xnsect=nsect
c
c  rotate ffts to new coordinate system
c  theta1 is angle of first component pos cw wrt geomagnetic north
c  theta2 is same for second component
c  phi is pos cw rotation angle of geomagnetic coordinate system
c  if original data are collected in orthogonal coordinate
c  system, set theta1=0, theta2=90
c
        if(nout.gt.1)then
c  first two ffts are assumed to be x and y and are rotated
c  if nout=3, third series is assumed to be z and is not rotated
          do 200 i=1,nfsect
            call rotate(nsect,theta1(1),theta2(1),phi(1),
     $                  ffts(1,1,i),ffts(1,2,i))
  200       continue
        endif
c
c  always rotate input ffts
c  if ninp=3, third series is assumed to be z and is not rotated
c
        if(ninp.gt.1)then
          do 202 i=1,nfsect
            call rotate(nsect,theta1(2),theta2(2),phi(2),
     $                  ffts(1,nout+1,i),ffts(1,nout+2,i))
  202       continue
        endif
c
c  always rotate remote reference ffts if nref.gt.1
c
        if(nref.gt.1)then
          do 204 j=1,nr3
            do 204 i=1,nfsect
              call rotate(nsect,theta1(2+j),theta2(2+j),phi(2+j),
     $                    ffts(1,nout+ninp+3*j-2,i),
     $                    ffts(1,nout+ninp+3*j-1,i))
  204         continue
          do 205 j=1,nr2
            do 205 i=1,nfsect
              call rotate(nsect,theta1(2+nr3+j),theta2(2+nr3+j),
     $                    phi(2+nr3+j),
     $                    ffts(1,nout+ninp+3*nr3+2*j-1,i),
     $                    ffts(1,nout+ninp+3*nr3+2*j,i))
  205         continue
        endif
c
c  compute degrees of freedom inflation factor
c
        dofscal=1.
        nm=max(0,min(nfft,(nptstot-nfft-npcs*abs(nar))/(nsect-1)))
        do 210 i=1,nsect-1
          sum=0.d0
          do 209 j=1,nfft
            if(j+i*nm.le.nfft)sum=sum+dpss(j)*dpss(j+i*nm)
  209       continue
  210     dofscal=dofscal+2.*(1-float(i/nsect))*sum**2
        dofscal=1./dofscal
c
c  begin second loop pivoting on frequency
c
        do 780 nf=1,nfsect
          frq=float(nf1+(nf-1)*nfinc)/xnfft
          do 215 k=1,nprej
            if(compare(frq,frej(k),4))goto 780
  215       continue
          freq(nfl)=frq
          write(19,*)' $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
          write(19,*)' period=',deltat/freq(nfl)
	  do 218 i=1,nsect
	      llo1(i)=.false.
	      lhi1(i)=.false.
              lwt1(i)=.true.
  218         wt1(i)=1.0
c
c  compute projected local magnetic field if required
c
          if(nrr.eq.0)goto 400
c
c  begin first stage
c
c coherence estimation by block averaging
c
          do 219 i=1,nsect/nblk
            c2flag(i)=.true.
	    do 219 j=1,ninp
  219         c2b(i,j)=1.0
          if(((nlev.ge.2).or.(c2threshb.gt.0.)).and.
     $       ((freq(nfl).ge.f2lo).and.(freq(nfl).le.f2hi)))then
            call bcoher2_b(nsect,ninp,nref,c2threshb,
     $                     ffts(1,nout+1,nf),c2flag,c2b,info)
            if(info.ne.0)then
              write(19,*)
     $          ' WARNING-singular matrix detected in',
     $          ' stage 1 coherence threshold at nsc=',nsc,' nf=',nf
            endif
          endif
          nsct=nsect
          do 221 i=1,nsect/nblk
            if(.not.c2flag(i))then
              do 220 ii=(i-1)*nblk+1,i*nblk
                nsct=nsct-1
  220           lwt1(ii)=.false.
            endif
  221       continue
          write(19,*)' stage 1 coherence threshold=',c2threshb
          if(c2threshb.gt.0.)
     $      write(19,*)nsect-nsct,' sections discarded'
          write(19,*)' initial number of data=',nsct
          if(nsct.eq.0)goto 780
          xnsct=nsct
c
c  begin third loop pivoting on input variable
c
	  do 385 nn=1,ninp
            write(19,*)' ++++++++++++++++++++++++++++++'
            write(19,*)' input series=',nn
c
c  compute section powers
c
	    fracthresh1(nn,nfl)=xnsct/xnsect
            do 224 i=1,nsect
              pwro(i)=2.*abs(ffts(i,nout+nn,nf))**2
              pwri(i)=0.
              do 223 j=1,nref
  223           pwri(i)=pwri(i)+abs(ffts(i,nout+ninp+j,nf))**2
	      pwri(i)=2.*pwri(i)/nref
  224         continue
c
c  compute nonrobust transfer function
c
            do 225 i=1,nsect
              wt(i)=0.
  225         if(lwt1(i))wt(i)=1.
            call transfer(nsect,ninp,nref,nn,wt,ffts(1,nout+1,nf),
     $                    b,aa,bb,info)
            if(info.eq.1)then
              write(19,*)' WARNING-singular matrix detected in',
     $                   ' nonrobust magnetic transfer section'
            endif
            write(19,*)
     $       ' initial transfer function indexed by reference variable'
            do 230 i=1,nref
              tbn(nn,i,nfl)=b(i)
  230         write(19,*)i,tbn(nn,i,nfl)
c
c  compute residuals
c
            call residual(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                    b,wt,nw,res,oresq)
c
c  compute jackknife estimate of transfer function standard error
c  compute hat matrix diagonal
c
            call hatdiag(nsect,nref,ffts(1,nout+ninp+1,nf),
     $                   wt,hat,info)
            if(info.ne.0)then
              write(19,*)
     $              ' WARNING-pivot block singular in nonrobust',
     $              ' magnetic transfer call to hatdiag info=',info
            endif
c
c  compute delete one solutions for transfer functions
c
            call delone1(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                   wt,aa,bb,hd1,info)
            if(info.ne.0)then
              write(19,*)
     $              ' WARNING-singular matrix detected in',
     $              ' nonrobust magnetic transfer jackknife'
            endif
c
c  form jackknife variance excluding estimates with zero weights
c
            do 240 ii=1,nref
              var=0.
              do 239 i=1,nsect
                if(wt(i).gt.0.)then
                  var=var+abs((1.-hat(i))*(tbn(nn,ii,nfl)-
     $                                     hd1(ii,i)))**2
                endif
  239           continue
  240         varn1(nn,ii,nfl)=xnsct*var/(xnsct-nref)
c
c  compute nonrobust complex coherence between measured and predicted output
c  variables
c
            call ccoher1(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                   tbn(nn,1,nfl),hd1,wt,gd1,gn1(nn,nfl),
     $                   gn1lo(nn,nfl),gn1hi(nn,nfl))
            dfn1(nn,nfl)=2.*dofscal*(nsct-ninp)
c
c output nonrobust diagnostics if requested
c
            if(abs(nlev).ge.1)then
              if(nfl.le.9)then
                cn(1:1)=char(48)
                cn(2:2)=char(48+nfl)
              elseif(nfl.le.99)then
                write(cn(1:2),'(i2)')nfl
              else
                write(cn,'(i3)')nfl
              endif
              call strlen(cn,m1,m2)
              filnam=ofil(l1:l2)//'.'//char(48+nn)//
     $               '.oq1.'//cn(m1:m2)
              call qqout(filnam,deltat/freq(nfl),nsect,nsct,nw,
     $                   nref,llo1,lhi1,wt,res,hat)
            endif
            if(abs(nlev).eq.3)then
              filnam=ofil(l1:l2)//'.'//char(48+nn)//'.or1.'//
     $               cn(m1:m2)
              call diagout(filnam,deltat/freq(nfl),nsect,nref,
     $                     res,pwro,pwri,wt,hat)
            endif
            if(uin.eq.-999.)then
              do 241 i=1,nref
  241           tbr(nn,i,nfl)=tbn(nn,i,nfl)
              goto 385
            endif
c
c begin robust transfer function computation
c
c
c  compute initial robust scale estimate using nonrobust residuals
c
            nr=0
	    do 245 i=1,nsect
	      xindex(i)=i
              if(wt(i).gt.0.)then
                nr=nr+1
                temp(nr)=res(i)
              endif
  245         continue
            if(nr.eq.0)goto 780
            call mad(nr,temp,scale)
c
c  huber weight regression section
c
c
c  compute number of hat matrix iterations
c
	    call ssort(hat,xindex,nsect,2)
	    khat=1
            if(ainuin.eq.-999.)then
              ainu=0.
              ainu1=0.
            elseif((ainuin.gt.0.).and.(ainuin.lt.1.))then
              ainu=nsct*qbeta(ainuin,nref,nsct-nref)/nref
	      ainu1=xnsct*hat(nsct)/nref
	      khat=max(log(ainu1/ainu)/log(2.)+1.,1.)
            elseif(.not.ilev)then
              ainu=nsct*qbeta(0.999,nref,nsct-nref)/nref
	      ainu1=xnsct*hat(nsct)/nref
	      khat=max(log(ainu1/ainu)/log(2.)+1.,1.)
            else
              ainu=abs(ainuin)
	      ainu1=xnsct*hat(nsct)/nref
	      khat=max(log(ainu1/ainu)/log(2.)+1.,1.)
            endif
            if(ainlin.eq.-999.)then
              ainl=0.
              ainl1=0.
            else
              do 256 k=1,nsect
  256           if(hat(k).ne.0.)goto 257
  257         ainl1=xnsct*hat(k)/nref
              if((ainlin.gt.0.).and.(ainlin.lt.1.))then
                ainl=nsct*qbeta(ainlin,nref,nsct-nref)/nref
              else
                ainl=abs(ainlin)
              endif
	      khat=max(float(khat),log(ainl/ainl1)/log(2.)+1.)
            endif
	    call ssort(xindex,hat,nsect,2)
	    write(19,*)' hat matrix iterations=',khat
            do 258 i=1,nsect
              what(i)=0.
  258         if(lwt1(i))what(i)=1.
            do 290 k=1,khat
c
c  begin hat matrix weight loop
c
              ainu1=max(ainu1/2.,ainu)
              ainl1=min(2.*ainl1,ainl)
              write(19,*)' lower cutoff=',ainl1,' upper cutoff=',ainu1
              do 280 niter=1,10
c
c begin huber regression weight loop
c
                write(19,*)' huber iteration=',niter,' scale=',scale
                if((ainlin.ne.-999.).or.(ainuin.ne.-999.))
     $            call hatwt(nsect,nref,ainl1,ainu1,llo1,lhi1,hat,what)
                call huberwt(nsect,scale,lwt1,res,w)
		do 260 i=1,nsect
  260             wt(i)=what(i)*w(i)
c
c  compute robust transfer function
c
                call transfer(nsect,ninp,nref,nn,wt,
     $                        ffts(1,nout+1,nf),b,aa,bb,info)
                if(info.eq.1)then
                  write(19,*)
     $                  ' WARNING-singular matrix detected in',
     $                  '  huber magnetic transfer section'
                  write(19,*)' niter=',niter
                endif
c
c  compute residuals
c
                call residual(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                        b,wt,nw,res,resq)
c
c  compute robust scale estimate
c
                nr=0
                do 276 i=1,nsect
                  if(wt(i).gt.0.)then
                    nr=nr+1
                    temp(nr)=res(i)
                  endif
  276             continue
                if(nr.eq.0)goto 780
                call mad(nr,temp,scale)
                if((ainlin.ne.-999.).or.(ainuin.ne.-999.))then
c
c  compute hat matrix diagonal
c
                  call hatdiag(nsect,nref,ffts(1,nout+ninp+1,nf),
     $                         wt,hat,info)
                  if(info.ne.0)then
                    write(19,*)
     $                ' WARNING-pivot block singular in huber magnetic',
     $                ' transfer call to hatdiag'
                    write(19,*)' info=',info,' niter=',niter
                  endif
                endif
                write(19,*)' normalized residual sum of squares=',
     $                       abs(resq-oresq)/abs(resq)
                if(abs(resq-oresq).lt.2.*eps*abs(resq))goto 285
  280           oresq=resq
              write(19,*)
     $              ' WARNING-huber iterations did not converge',
     $              ' at hat matrix iteration=',k
  285         oresq=resq
              write(19,*)' number of data=',nw
              write(19,*)
     $         ' huber transfer function indexed by reference variable'
              do 290 i=1,nref
                ob(i)=b(i)
  290           write(19,*)i,cmplx(b(i))
c
c  thomson weight regression section
c
            if(.not.ilev)then
              if((uin.eq.0.).or.(abs(uin).ge.1.))then
                u0=sqrt(2.*log(xnsct/(abs(uin)+0.5)))
              elseif((abs(uin).gt.0.).and.(abs(uin).lt.1.))then
                u0=sqrt(-2.*log(1.-abs(uin)))
              endif
            else
              if((uin.eq.0.).or.(uin.ge.1.))then
                u0=sqrt(2.*log(xnsct/(uin+0.5)))
              elseif((uin.gt.0.).and.(uin.lt.1.))then
                u0=sqrt(-2.*log(1.-uin))
              else
                u0=abs(uin)
              endif
            endif
            if((ainlin.ne.-999.).or.(ainuin.ne.-999.))
     $        call hatwt(nsect,nref,ainl,ainu,llo1,lhi1,hat,what)
            call expwt(nsect,u0,scale,lwt1,res,w)
            do 300 i=1,nsect
  300         wt(i)=what(i)*w(i)
            do 340 niter=1,15
c
c  begin iteration on thomson regression weights
c
              write(19,*)
     $              ' thomson iteration=',niter,' scale=',scale
c
c  compute robust transfer function
c
              call transfer(nsect,ninp,nref,nn,wt,
     $                      ffts(1,nout+1,nf),b,aa,bb,info)
              if(info.eq.1)then
                write(19,*)
     $                ' WARNING-singular matrix detected in thomson',
     $                '  magnetic transfer section'
              endif
c
c  compute residuals
c
              call residual(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                      b,wt,nw,res,resq)
              write(19,*)' number of data=',nw
              write(19,*)
     $         ' thomson solution vector indexed by reference variable'
              do 315 i=1,nref
  315           write(19,*)i,cmplx(b(i))
              write(19,*)' normalized residual sum of squares=',
     $                     abs(resq-oresq)/abs(resq)
c
c  compute robust scale estimate
c
c              if(niter.le.2)then
              if(niter.eq.1)then
                nr=0
                do 320 i=1,nsect
                  if(wt(i).gt.0.)then
                    nr=nr+1
                    temp(nr)=res(i)
                  endif
  320             continue
                if(nr.eq.0)goto 780
                call mad(nr,temp,scale)
              endif
c
c  compute hat matrix diagonal
c
              call hatdiag(nsect,nref,ffts(1,nout+ninp+1,nf),
     $                     wt,hat,info)
              if(info.ne.0)then
                write(19,*)
     $                  ' WARNING-pivot block singular in',
     $                  ' thomson magnetic transfer call to hatdiag'
                write(19,*)' info=',info,' niter=',niter
              endif
              if((ainuin.ne.-999.).or.(ainlin.ne.-999.))
     $          call hatwt(nsect,nref,ainl,ainu,llo1,lhi1,hat,what)
	      call expwt(nsect,u0,scale,lwt1,res,w)
	      do 325 i=1,nsect
  325           wt(i)=what(i)*w(i)
c
c  compute jackknife estimate of transfer function standard error
c
c
c  compute delete one solutions for transfer functions
c
              call delone1(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                     wt,aa,bb,hd1,info)
	      if(info.ne.0)then
	        write(19,*)
     $              ' WARNING-singular matrix detected in stage 1',
     $              ' robust jackknife niter=',niter
              endif
c
c  form jackknife variance excluding values with zero weights
c
              do 330 ii=1,nref
	        var=0.
	        do 327 i=1,nsect
	          if(wt(i).gt.0.)then
		    var=var+abs((1.-hat(i))*
     $                      (b(ii)-hd1(ii,i)))**2
                  endif
  327         continue
  330         varr1(nn,ii,nfl)=nw*var/(nw-ninp)
	      cflag1=.true.
	      do 335 i=1,nref
		sdev=sqrt(varr1(nn,i,nfl))
		if((abs(dreal(ob(i)-b(i))).ge.eps1*sdev).or.
     $             (abs(dimag(ob(i)-b(i))).ge.eps1*sdev))
     $             cflag1=.false.
  335           ob(i)=b(i)
              cflag2=.false.
	      if((abs(real(resq-oresq)).lt.eps*real(resq)).and.
     $           (abs(aimag(resq-oresq)).lt.eps*aimag(resq)))
     $           cflag2=.true.
	      if(cflag1.or.cflag2)goto 345
  340         oresq=resq
	    convflag1(nn,nfl)=1
	    write(19,*)' WARNING-thomson iterations did not converge'
  345       if(cflag1)write(19,*)' convergence on response'
	    if(cflag2)write(19,*)' convergence on residuals'
  350       do 360 i=1,nref
  360         tbr(nn,i,nfl)=b(i)
c
c  compute complex coherence between measured and predicted output variables
c
            call ccoher1(nsect,ninp,nref,nn,ffts(1,nout+1,nf),
     $                   tbr(nn,1,nfl),hd1,wt,gd1,gr1(nn,nfl),
     $                   gr1lo(nn,nfl),gr1hi(nn,nfl))
	    wtsum=0.
	    do 370 i=1,nsect
  370         wtsum=wtsum+wt(i)
	    dfr1(nn,nfl)=2.*dofscal*(wtsum-ninp)
	    fracdat1(nn,nfl)=fracthresh1(nn,nfl)*wtsum/xnsct
c
c  output robust diagnostics as required
c
	    if(abs(nlev).ge.1)then
	      filnam=ofil(l1:l2)//'.'//char(48+nn)//
     $               '.fq1.'//cn(m1:m2)
              call hatdiag(nsect,nref,ffts(1,nout+ninp+1,nf),
     $                     wt,hat,info)
              if(info.ne.0)then
                write(19,*)
     $                ' WARNING-pivot block singular in',
     $                ' thomson magnetic transfer call to hatdiag'
                write(19,*)' info=',info,' niter=',niter
              endif
              call qqout(filnam,deltat/freq(nfl),nsect,nsct,nw,
     $                   nref,llo1,lhi1,wt,res,hat)
	    endif
	    if(abs(nlev).eq.3)then
	      filnam=ofil(l1:l2)//'.'//char(48+nn)//'.fr1.'//
     $               cn(m1:m2)
	      call diagout(filnam,deltat/freq(nfl),nsect,nref,
     $                     res,pwro,pwri,wt,hat)
            endif
 	    do 381 i=1,nsect
  381         wt1(i)=min(wt1(i),wt(i))
  385       continue
c
c  compute projected local magnetic field
c
        do 395 k=1,nsect
	      do 390 i=1,ninp
	        b(i)=0.
            do 390 j=1,nref
  390         b(i)=b(i)+ffts(k,nout+ninp+j,nf)*tbr(i,j,nfl)
            do 395 i=1,ninp
              ffts(k,nout+ninp+i,nf)=b(i)
  395         ffts(k,nout+i,nf)=ffts(k,nout+ninp+i,nf)
  400     continue
c
c  begin second stage
c
c coherence estimation by block averaging
c
          do 405 i=1,nsect/nblk
            c2flag(i)=.true.
            c2flag1(i)=.true.
            do 405 j=1,nout
  405         c2e(i,j)=1.0
          if(((nlev.ge.2).or.(c2threshe.gt.0.).or.
     $       (c2threshe1.gt.0.)).and.
     $       ((freq(nfl).ge.f2lo).and.(freq(nfl).le.f2hi)))then
            n1=min(nout,2)
            if((nout.eq.3).and.(nz.gt.0))n1=3
            call bcoher2_e(nsect,n1,nout,ninp,c2threshe,
     $                     ffts(1,1,nf),c2flag,c2e,info)
            if(info.ne.0)then
              write(19,*)
     $          ' WARNING-singular matrix detected in',
     $          ' stage 2 coherence threshold at nsc=',nsc,' nf=',nf
            endif
            if(nz.gt.0)then
              do 406 i=1,nsect/nblk
  406           c2flag1(i)=c2flag(i)
            endif
            if((nout.eq.3).and.(nz.eq.0))then
              call bcoher2_e(nsect,1,1,ninp,c2threshe1,
     $                       ffts(1,nout,nf),c2flag1,c2e(1,nout),
     $                       info)
              if(info.ne.0)then
                write(19,*)
     $            ' WARNING-singular matrix detected in',
     $            ' stage 2 coherence threshold at nsc=',nsc,' nf=',nf
              endif
            endif
          endif
          write(19,*)' stage 2 coherence threshold=',c2threshe
          if(nout.eq.3)
     $       write(19,*)' stage 2 3rd output coherence threshold=',
     $                    c2threshe1
c
c  begin third loop pivoting on output variable
c
          do 770 nn=1,nout
            if(nfl.le.9)then
              cn(1:1)=char(48)
              cn(2:2)=char(48+nfl)
            elseif(nfl.le.99)then
              write(cn(1:2),'(i2)')nfl
            else
              write(cn,'(i3)')nfl
            endif
            write(19,*)' #####################################'
            write(19,*)' output series=',nn
            if((nn.lt.3).or.(nz.gt.0))then
              do 410 i=1,nsect
  410           lwt(i)=lwt1(i)
              do 412 i=1,nsect/nblk
                if(.not.c2flag(i))then
                  do 411 ii=(i-1)*nblk+1,i*nblk
  411               lwt(ii)=.false.
                endif
  412           continue
            else
              do 415 i=1,nsect
  415           lwt(i)=lwt1(i)
              do 417 i=1,nsect/nblk
                if(.not.c2flag1(i))then
                  do 416 ii=(i-1)*nblk+1,i*nblk
  416               lwt(ii)=.false.
                endif
  417           continue
            endif
            nsct=0
            do 418 i=1,nsect
  418         if((wt1(i).gt.0.).and.lwt(i))nsct=nsct+1
            xnsct=nsct
            if((c2threshe.gt.0.).or.(c2threshe1.gt.0.).or.(nrr.eq.1))
     $        write(19,*)nsect-nsct,' sections discarded'
            write(19,*)' initial number of data=',nsct
	    fracthresh2(nn,nfl)=xnsct/xnsect
            if(nsct.eq.0)goto 770
c
c  compute section powers
c
            do 420 i=1,nsect
              pwro(i)=2.*abs(ffts(i,nn,nf))**2
              pwri(i)=0.
              do 419 j=1,ninp
  419           pwri(i)=pwri(i)+abs(ffts(i,nout+j,nf))**2
	      pwri(i)=2.*pwri(i)/ninp
  420         continue
c
c  compute nonrobust transfer function
c
            do 421 i=1,nsect
              wt(i)=wt1(i)
  421       if(.not.lwt(i))wt(i)=0.
            call response(nsect,nout,ninp,nn,wt,ffts(1,1,nf),b,aa,
     $                    bb,info)
            if(info.eq.1)then
              write(19,*)' WARNING-singular matrix detected in',
     $                     ' l2 nonrobust response section '
            endif
            write(19,*)
     $            ' initial solution vector indexed by input variable'
            do 422 i=1,ninp
              hn(nn,i,nfl)=b(i)
  422         write(19,*)i,hn(nn,i,nfl)
c
c  compute residuals
c
            call residual(nsect,nout,ninp,nn,ffts(1,1,nf),
     $                    b,wt,nw,res,oresq)
c
c  compute jackknife estimate of transfer function standard error
c  compute hat matrix diagonal
c
            call hatdiag(nsect,ninp,ffts(1,nout+ninp+1,nf),
     $                   wt,hat,info)
            if(info.ne.0)then
              write(19,*)
     $              ' WARNING-pivot block singular in nonrobust',
     $              ' response call to hatdiag  info=',info
            endif
c
c  compute delete one solutions for transfer functions
c
            call delone2(nsect,nout,ninp,nn,ffts(1,1,nf),wt,aa,
     $                bb,hd1,info)
            if(info.ne.0)then
              write(19,*)
     $              ' WARNING-singular matrix detected in',
     $              ' nonrobust jackknife'
            endif
c
c  form jackknife variance excluding estimates with zero weights
c
            do 430 ii=1,ninp
              var=0.
              do 429 i=1,nsect
                if(wt(i).gt.0.)then
                  var=var+abs((1.-hat(i))*(hn(nn,ii,nfl)-
     $                                     hd1(ii,i)))**2
                endif
  429           continue
  430         varn2(nn,ii,nfl)=xnsct*var/(xnsct-ninp)
c
c  compute nonrobust complex coherence between measured and predicted output
c  variables
c
            call ccoher2(nsect,nout,ninp,nn,ffts(1,1,nf),hn(nn,1,nfl),
     $                  hd1,wt,gd1,gn2(nn,nfl),gn2lo(nn,nfl),
     $                  gn2hi(nn,nfl))
            dfn2(nn,nfl)=2.*dofscal*(nsct-ninp)
c
c output nonrobust diagnostics if requested
c
            if(abs(nlev).ge.1)then
              if(nfl.le.9)then
                cn(1:1)=char(48)
                cn(2:2)=char(48+nfl)
              elseif(nfl.le.99)then
                write(cn(1:2),'(i2)')nfl
              else
                write(cn,'(i3)')nfl
              endif
              call strlen(cn,m1,m2)
              filnam=ofil(l1:l2)//'.'//char(48+nn)//
     $               '.oq2.'//cn(m1:m2)
              call qqout(filnam,deltat/freq(nfl),nsect,nsct,nw,
     $                   ninp,llo1,lhi1,wt,res,hat)
            endif
            if(abs(nlev).eq.3)then
              filnam=ofil(l1:l2)//'.'//char(48+nn)//'.or2.'//
     $               cn(m1:m2)
              call diagout(filnam,deltat/freq(nfl),nsect,ninp,
     $                     res,pwro,pwri,wt,hat)
            endif
c
c begin robust transfer function computation
c
c
c  compute initial robust scale estimate using nonrobust residuals
c
            if(uin.eq.-999.)goto 760
            nr=0
            do 490 i=1,nsect
              xindex(i)=i
              if(wt(i).gt.0.)then
                nr=nr+1
                temp(nr)=res(i)
              endif
  490         continue
            if(nr.eq.0)goto 770
            call mad(nr,temp,scale)
c
c  huber weight regression section
c
c
c  compute number of hat matrix iterations
c
            call ssort(hat,xindex,nsect,2)
            do 495 k=1,nsect
  495         if(hat(k).ne.0.)goto 496
  496       khat=1
            if(ainuin.eq.-999.)then
              ainu=0.
              ainu1=0.
            elseif((ainuin.gt.0.).and.(ainuin.lt.1.))then
              ainu=nsct*qbeta(ainuin,ninp,nsct-ninp)/ninp
	      ainu1=xnsct*hat(nsct+k-1)/ninp
	      khat=max(log(ainu1/ainu)/log(2.)+1.,1.)
            else
              ainu=abs(ainuin)
	      ainu1=xnsct*hat(nsct+k-1)/ninp
	      khat=max(log(ainu1/ainu)/log(2.)+1.,1.)
            endif
            if(ainlin.eq.-999.)then
              ainl=0.
              ainl1=0.
            else
              ainl1=xnsct*hat(k)/ninp
              if((ainlin.gt.0.).and.(ainlin.lt.1.))then
                ainl=nsct*qbeta(ainlin,ninp,nsct-ninp)/ninp
              else
                ainl=abs(ainlin)
              endif
	      khat=max(float(khat),log(ainl/ainl1)/log(2.)+1.)
            endif
            call ssort(xindex,hat,nsect,2)
            write(19,*)' hat matrix iterations=',khat
            do 500 i=1,nsect
              llo(i)=llo1(i)
              lhi(i)=lhi1(i)
              what(i)=0.
  500         if((wt1(i).gt.0.).and.lwt(i))what(i)=1.
            do 560 k=1,khat
c
c  begin hat matrix weight loop
c
              ainu1=max(ainu1/2.,ainu)
              ainl1=min(2.*ainl1,ainl)
              write(19,*)' lower cutoff=',ainl1,' upper cutoff=',ainu1
              do 550 niter=1,10
c
c begin huber regression weight loop
c
                write(19,*)' huber iteration=',niter,' scale=',scale
                if((ainlin.ne.-999.).or.(ainuin.ne.-999.))
     $            call hatwt(nsect,ninp,ainl1,ainu1,llo,lhi,hat,what)
                call huberwt(nsect,scale,lwt,res,w)
                do 510 i=1,nsect
  510             wt(i)=wt1(i)*what(i)*w(i)
c
c  compute robust transfer function
c
                call response(nsect,nout,ninp,nn,wt,ffts(1,1,nf),b,
     $                        aa,bb,info)
                if(info.eq.1)then
                  write(19,*)
     $                  ' WARNING-singular matrix detected in',
     $                  '  huber section'
                  write(19,*)' niter=',niter
                endif
c
c  compute residuals
c
                call residual(nsect,nout,ninp,nn,ffts(1,1,nf),
     $                        b,wt,nw,res,resq)
c
c  compute robust scale estimate
c
                nr=0
                do 540 i=1,nsect
                  if(wt(i).gt.0.)then
                    nr=nr+1
                    temp(nr)=res(i)
                  endif
  540             continue
                if(nr.eq.0)goto 770
                call mad(nr,temp,scale)
                if((ainlin.ne.-999.).or.(ainuin.ne.-999.))then
c
c  compute hat matrix diagonal
c
                  call hatdiag(nsect,ninp,ffts(1,nout+ninp+1,nf),
     $                         wt,hat,info)
                  if(info.ne.0)then
                    write(19,*)
     $                ' WARNING-pivot block singular in huber call',
     $                ' to hatdiag'
                    write(19,*)' info=',info,' niter=',niter
                  endif
                endif
                write(19,*)' normalized residual sum of squares=',
     $                       abs(resq-oresq)/abs(resq)
                if(abs(resq-oresq).lt.2.*eps*abs(resq))
     $             goto 555
  550           oresq=resq
              write(19,*)
     $              ' WARNING-huber iterations did not converge',
     $              ' at hat matrix iteration=',k
  555         oresq=resq
              write(19,*)' number of data=',nw
              write(19,*)
     $              ' huber solution vector indexed by input variable'
              do 560 i=1,ninp
                ob(i)=b(i)
  560           write(19,*)i,cmplx(b(i))
c
c  thomson weight regression section
c
            if(.not.ilev)then
              if((uin.eq.0.).or.(abs(uin).ge.1.))then
                u0=sqrt(2.*log(xnsct/(abs(uin)+0.5)))
              elseif((abs(uin).gt.0.).and.(abs(uin).lt.1.))then
                u0=sqrt(-2.*log(1.-abs(uin)))
              endif
            else
              if((uin.eq.0.).or.(uin.ge.1.))then
                u0=sqrt(2.*log(xnsct/(uin+0.5)))
              elseif((uin.gt.0.).and.(uin.lt.1.))then
                u0=sqrt(-2.*log(1.-uin))
              else
                u0=abs(uin)
              endif
            endif
            if((ainlin.ne.-999.).or.(ainuin.ne.-999.))
     $        call hatwt(nsect,ninp,ainl,ainu,llo,lhi,hat,what)
            call expwt(nsect,u0,scale,lwt,res,w)
            do 565 i=1,nsect
  565         wt(i)=wt1(i)*what(i)*w(i)
            do 600 niter=1,15
c
c  begin iteration on thomson regression weights
c
              write(19,*)
     $              ' thomson iteration=',niter,' scale=',scale
c
c  compute robust transfer function
c
              call response(nsect,nout,ninp,nn,wt,ffts(1,1,nf),b,
     $                      aa,bb,info)
              if(info.eq.1)then
                write(19,*)
     $                ' WARNING-singular matrix detected in thomson',
     $                '  section'
              endif
c
c  compute residuals
c
              call residual(nsect,nout,ninp,nn,ffts(1,1,nf),
     $                      b,wt,nw,res,resq)
              write(19,*)' number of data=',nw
              write(19,*)
     $              ' thomson solution vector indexed by input variable'
              do 570 i=1,ninp
  570           write(19,*)i,cmplx(b(i))
              write(19,*)' normalized residual sum of squares=',
     $                     abs(resq-oresq)/abs(resq)
c
c  compute robust scale estimate
c
c              if(niter.le.2)then
              if(niter.eq.1)then
                nr=0
                do 575 i=1,nsect
                  if(wt(i).gt.0.)then
                    nr=nr+1
                    temp(nr)=res(i)
                  endif
  575             continue
                if(nr.eq.0)goto 770
                call mad(nr,temp,scale)
              endif
c
c  compute hat matrix diagonal
c
              call hatdiag(nsect,ninp,ffts(1,nout+ninp+1,nf),
     $                     wt,hat,info)
              if(info.ne.0)then
                write(19,*)
     $                ' WARNING-pivot block singular in',
     $                ' thomson call to hatdiag'
                write(19,*)' info=',info,' niter=',niter
              endif
              if((ainuin.ne.-999.).or.(ainlin.ne.-999.))
     $          call hatwt(nsect,ninp,ainl,ainu,llo,lhi,hat,what)
	      call expwt(nsect,u0,scale,lwt,res,w)
	      do 580 i=1,nsect
  580           wt(i)=wt1(i)*what(i)*w(i)
c
c  compute jackknife estimate of transfer function standard error
c
c
c  compute delete one solutions for transfer functions
c
              call delone2(nsect,nout,ninp,nn,ffts(1,1,nf),wt,aa,
     $                     bb,hd1,info)
	      if(info.ne.0)then
	        write(19,*)
     $              ' WARNING-singular matrix detected in robust',
     $              ' jackknife niter=',niter
              endif
c
c  form jackknife variance excluding values with zero weights
c
              do 590 ii=1,ninp
	        var=0.
	        nw=0
	        do 585 i=1,nsect
	          if(wt(i).gt.0.)then
	            nw=nw+1
		    var=var+abs((1.-hat(i))*(b(ii)-hd1(ii,i)))**2
                  endif
  585         continue
  590         varr2(nn,ii,nfl)=nw*var/(nw-ninp)
	      cflag1=.true.
	      do 595 i=1,ninp
		sdev=sqrt(varr2(nn,i,nfl))
		if((abs(dreal(ob(i)-b(i))).ge.eps1*sdev).or.
     $             (abs(dimag(ob(i)-b(i))).ge.eps1*sdev))
     $             cflag1=.false.
  595           ob(i)=b(i)
		cflag2=.false.
		if((abs(real(resq-oresq)).lt.eps*real(resq)).and.
     $             (abs(aimag(resq-oresq)).lt.eps*aimag(resq)))
     $             cflag2=.true.
c	      if((niter.gt.1).and.(cflag1.or.cflag2))goto 610
	      if(cflag1.or.cflag2)goto 610
  600         oresq=resq
	    convflag2(nn,nfl)=1
	    write(19,*)' WARNING-thomson iterations did not converge'
  610       if(cflag1)write(19,*)' convergence on response'
	    if(cflag2)write(19,*)' convergence on residuals'
	    do 620 i=1,ninp
  620         hr(nn,i,nfl)=b(i)
c
c  compute complex coherence between measured and predicted output variables
c
            call ccoher2(nsect,nout,ninp,nn,ffts(1,1,nf),hr(nn,1,nfl),
     $                  hd1,wt,gd1,gr2(nn,nfl),gr2lo(nn,nfl),
     $                  gr2hi(nn,nfl))
	    wtsum=0.
	    do 640 i=1,nsect
  640         wtsum=wtsum+wt(i)
	    dfr2(nn,nfl)=2.*dofscal*(wtsum-ninp)
	    fracdat2(nn,nfl)=fracthresh2(nn,nfl)*wtsum/xnsct
	    if((fracdat2(nn,nfl).ne.0.).and.(nlev.lt.0))then
	      write(32+nn)freq(nfl),nsect
	      write(32+nn)(hat(i),wt(i),(hd1(j,i),j=1,ninp),
     $                     i=1,nsect)
	    endif
c
c  output robust diagnostics as required
c
	    if(abs(nlev).ge.1)then
	      filnam=ofil(l1:l2)//'.'//char(48+nn)//
     $               '.fq2.'//cn(m1:m2)
              call qqout(filnam,deltat/freq(nfl),nsect,nsct,nw,
     $                   ninp,llo,lhi,wt,res,hat)
	    endif
	    if(abs(nlev).eq.3)then
	      filnam=ofil(l1:l2)//'.'//char(48+nn)//'.fr2.'//
     $               cn(m1:m2)
              call diagout(filnam,deltat/freq(nfl),nsect,ninp,
     $                     res,pwro,pwri,wt,hat)
	    endif
  760       if(abs(nlev).ge.2)then
	      filnam=ofil(l1:l2)//'.'//char(48+nn)//
     $               '.wt.'//cn(m1:m2)
	      open (unit=40,file=filnam,status='unknown')
	      if(nrr.eq.0)then
		do 764 ii=1,nsect/nblk
		  do 762 i=1,nblk
                    k=(ii-1)*nblk+i
		    write(40,999)deltat*(xnfft/2.+float(k-1)*noff),
     $                           c2e(ii,nn),w(k),what(k)
  762               continue
  764             continue
	      else
		do 768 ii=1,nsect/nblk
		  do 766 i=1,nblk
                    k=(ii-1)*nblk+i
	            write(40,999)deltat*(xnfft/2.+float(k-1)*noff),
     $                           c2e(ii,nn),(c2b(ii,j),j=1,ninp),
     $                           wt1(k),wt(k)
  766               continue
  768             continue
	      endif
	      close (unit=40)
	    endif
  770       continue
	  nfl=nfl+1
  780     continue
	nfft=nfft/nsctinc
  790   continue
c
c  output section
c
      do 902 n=1,nout
        if(nlev.lt.0)close(unit=32+n)
  902   continue
c
c sort output frequencies into ascending order
c
      do 905 i=1,nfl-1
  905   xindex(i)=i
      call ssort(freq,xindex,nfl-1,2)
      do 906 i=1,nfl-1
  906   nfr(i)=xindex(i)
      do 950 n=1,nout
	do 930 j=1,ninp
	  postfix='.'//char(48+n)//'n'//char(48+j)//'.rf'
	  call strlen(postfix,m1,m2)
	  filnam=ofil(l1:l2)//postfix(m1:m2)
	  open (unit=10,file=filnam,status='unknown')
	  do 910 i=1,nfl-1
	    ii=nfr(i)
	    if(fracthresh2(n,ii).ne.0.)then
	      write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                     real(hn(n,j,ii)),aimag(hn(n,j,ii)),
     $                     sqrt(varn2(n,j,ii)),dfn2(n,ii),
     $                     fracthresh2(n,ii)
	    endif
  910       continue
	  close (unit=10)
	  if(uin.ne.-999.)then
	    postfix='.'//char(48+n)//'r'//char(48+j)//'.rf'
	    filnam=ofil(l1:l2)//postfix(m1:m2)
            open (unit=10,file=filnam,status='unknown')
	    do 915 i=1,nfl-1
	      ii=nfr(i)
	      if(fracdat2(n,ii).ne.0.)then
		write(10,1000)deltat/freq(ii),freq(ii)/deltat,
     $                        real(hr(n,j,ii)),aimag(hr(n,j,ii)),
     $                        sqrt(varr2(n,j,ii)),dfr2(n,ii),
     $                        fracthresh2(n,ii),fracdat2(n,ii),
     $                        convflag2(n,ii)
	      endif
  915         continue
	    close (unit=10)
	  endif
	  if(n.le.2)then
	    postfix='.'//char(48+n)//'n'//char(48+j)//'.rp'
	    filnam=ofil(l1:l2)//postfix(m1:m2)
	    open (unit=10,file=filnam,status='unknown')
	    do 920 i=1,nfl-1
	      ii=nfr(i)
	      if(fracthresh2(n,ii).ne.0.)then
                sdev=sqrt(varn2(n,j,ii))
                if(2.24*sdev/abs(hn(n,j,ii)).lt.1.)then
                  write(10,1000)deltat/freq(ii),freq(ii)/deltat,
     $                   0.2*deltat/freq(ii)*abs(hn(n,j,ii))**2,
     $                   0.896*deltat/freq(ii)*abs(hn(n,j,ii))*
     $                   sdev,atan2(aimag(hn(n,j,ii)),
     $                   real(hn(n,j,ii)))/rad,
     $                   asin(2.24*sdev/abs(hn(n,j,ii)))/rad,
     $                   dfn2(n,ii),fracthresh2(n,ii)
                else
                  write(10,1000)deltat/freq(ii),freq(ii)/deltat,
     $                   0.2*deltat/freq(ii)*abs(hn(n,j,ii))**2,
     $                   0.896*deltat/freq(ii)*abs(hn(n,j,ii))*
     $                   sdev,atan2(aimag(hn(n,j,ii)),
     $                   real(hn(n,j,ii)))/rad,180.0,
     $                   dfn2(n,ii),fracthresh2(n,ii)
                endif
	      endif
  920         continue
	    close (unit=10)
	  endif
	  if((uin.ne.-999.).and.(n.le.2))then
	    postfix='.'//char(48+n)//'r'//char(48+j)//'.rp'
	    filnam=ofil(l1:l2)//postfix(m1:m2)
	    open (unit=10,file=filnam,status='unknown')
	    do 925 i=1,nfl-1
	      ii=nfr(i)
	      if(fracdat2(n,ii).ne.0.)then
                sdev=sqrt(varr2(n,j,ii))
                if(2.24*sdev/abs(hr(n,j,ii)).lt.1.)then
                  write(10,1001)deltat/freq(ii),freq(ii)/deltat,
     $                     0.2*deltat/freq(ii)*abs(hr(n,j,ii))**2,
     $                     0.896*deltat/freq(ii)*abs(hr(n,j,ii))*
     $                     sdev,atan2(aimag(hr(n,j,ii)),
     $                     real(hr(n,j,ii)))/rad,
     $                     asin(2.24*sdev/abs(hr(n,j,ii)))/rad,
     $                     dfr2(n,ii),fracthresh2(n,ii),
     $                     fracdat2(n,ii),convflag2(n,ii)
                else
                  write(10,1001)deltat/freq(ii),freq(ii)/deltat,
     $                     0.2*deltat/freq(ii)*abs(hr(n,j,ii))**2,
     $                     0.896*deltat/freq(ii)*abs(hr(n,j,ii))*
     $                     sdev,atan2(aimag(hr(n,j,ii)),
     $                     real(hr(n,j,ii)))/rad,
     $                     180.,dfr2(n,ii),fracthresh2(n,ii),
     $                     fracdat2(n,ii),convflag2(n,ii)
                endif
	      endif
  925         continue
	      close (unit=10)
	    endif
  930     continue
	  postfix='.'//char(48+n)//'n.2c2'
	  call strlen(postfix,m1,m2)
	  filnam=ofil(l1:l2)//postfix(m1:m2)
	  open (unit=10,file=filnam,status='unknown')
	  do 935 i=1,nfl-1
	    ii=nfr(i)
	    if(fracthresh2(n,ii).ne.0.)then
	      write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                   abs(gn2(n,ii))**2,phase(gn2(n,ii)),
     $                   1.-.05**(2./(dfn2(n,ii)-2.)),
     $                   gn2lo(n,ii)**2,gn2hi(n,ii)**2,dfn2(n,ii),
     $                   fracthresh2(n,ii)
	    endif
  935       continue
	  close (unit=10)
	  if(uin.ne.-999.)then
            postfix='.'//char(48+n)//'r.2c2'
	    filnam=ofil(l1:l2)//postfix(m1:m2)
	    open (unit=10,file=filnam,status='unknown')
	    do 940 i=1,nfl-1
	      ii=nfr(i)
	      if(fracdat2(n,ii).ne.0.)then
	        write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                     abs(gr2(n,ii))**2,phase(gr2(n,ii)),
     $                     1.-.05**(2./(dfr2(n,ii)-2.)),
     $                     gr2lo(n,ii)**2,gr2hi(n,ii)**2,dfr2(n,ii),
     $                     fracthresh2(n,ii),fracdat2(n,ii)
	      endif
  940         continue
	    close (unit=10)
	  endif
  950   continue
      if(nrr.eq.1)then
	do 980 n=1,ninp
	  do 955 j=1,nref
	    postfix='.'//char(48+n)//'n'//char(48+j)//'.tf'
	    call strlen(postfix,m1,m2)
	    filnam=ofil(l1:l2)//postfix(m1:m2)
	    open (unit=10,file=filnam,status='unknown')
	    do 955 i=1,nfl-1
	      ii=nfr(i)
	      if(fracthresh1(n,ii).ne.0.)then
	        write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                       real(tbn(n,j,ii)),aimag(tbn(n,j,ii)),
     $                       sqrt(varn1(n,j,ii)),dfn1(n,ii),
     $                       fracthresh1(n,ii)
	      endif
  955         continue
	    close (unit=10)
	    if(uin.ne.-999.)then
	      do 960 j=1,nref
	        postfix='.'//char(48+n)//'r'//char(48+j)//'.tf'
	        filnam=ofil(l1:l2)//postfix(m1:m2)
                open (unit=10,file=filnam,status='unknown')
	        do 960 i=1,nfl-1
	          ii=nfr(i)
	          if(fracdat1(n,ii).ne.0.)then
		    write(10,1000)deltat/freq(ii),freq(ii)/deltat,
     $                            real(tbr(n,j,ii)),aimag(tbr(n,j,ii)),
     $                            sqrt(varr1(n,j,ii)),dfr1(n,ii),
     $                            fracthresh1(n,ii),fracdat1(n,ii),
     $                            convflag1(n,ii)
	          endif
  960             continue
	      close (unit=10)
	    endif
	    postfix='.'//char(48+n)//'n.1c2'
	    call strlen(postfix,m1,m2)
	    filnam=ofil(l1:l2)//postfix(m1:m2)
	    open (unit=10,file=filnam,status='unknown')
	    do 965 i=1,nfl-1
	      ii=nfr(i)
	      if(fracthresh1(n,ii).ne.0.)then
	        write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                       abs(gn1(n,ii))**2,phase(gn1(n,ii)),
     $                       1.-.05**(2./(dfn1(n,ii)-2.)),
     $                       gn1lo(n,ii)**2,gn1hi(n,ii)**2,dfn1(n,ii),
     $                       fracthresh1(n,ii)
	      endif
  965         continue
	    close (unit=10)
	    if(uin.ne.-999.)then
              postfix='.'//char(48+n)//'r.1c2'
	      filnam=ofil(l1:l2)//postfix(m1:m2)
	      open (unit=10,file=filnam,status='unknown')
	      do 970 i=1,nfl-1
	        ii=nfr(i)
	        if(fracdat1(n,ii).ne.0.)then
	          write(10,999)deltat/freq(ii),freq(ii)/deltat,
     $                         abs(gr1(n,ii))**2,phase(gr1(n,ii)),
     $                         1.-.05**(2./(dfr1(n,ii)-2.)),
     $                         gr1lo(n,ii)**2,gr1hi(n,ii)**2,dfr1(n,ii),
     $                         fracthresh1(n,ii),fracdat1(n,ii)
	        endif
  970           continue
	      close (unit=10)
            endif
  980     continue
        endif
        if(uin.ne.-999)then
          write(18,'(2a)')'>LATITUDE  = ',cvar(1)
          write(18,'(2a)')'>LONGITUDE = ',cvar(2)
          write(18,'(2a)')'>ELEVATION = ',cvar(3)
          write(18,'(a,g15.7)')'>AZIMUTH   = ',phi(1)
          write(18,'(a)')ofil
          do 985 n=1,nout
	    do 985 j=1,2
              write(18,'(a)')jlabel2(n,j)
	      write(18,'(i2)')nfl-1
              do 985 i=nfl-1,1,-1
                ii=nfr(i)
	        if(fracdat2(n,ii).ne.0.)then
		  write(18,999)deltat/freq(ii),real(hr(n,j,ii)),
     $			       aimag(hr(n,j,ii)),sqrt(varr2(n,j,ii)),
     $                         fracthresh2(n,ii),fracdat2(n,ii)
	        else
	          write(18,999)-999.,-999.,-999.,
     $                         -999.,-999.,-999.
	        endif
 985	        continue
          do 987 n=1,2
    	    do 987 j=1,2
	      write(18,'(a)')jlabel1(n,j)
	      write(18,'(i2)')nfl-1
              do 987 i=nfl-1,1,-1
                ii=nfr(i)
                if(fracdat2(n,ii).ne.0.)then
		  rho=0.2*deltat/freq(ii)*abs(hr(n,j,ii))**2
                  ph=atan2(aimag(hr(n,j,ii)),real(hr(n,j,ii)))/rad
                  sdev=sqrt(varr2(n,j,ii))
		  drho=0.4*deltat/freq(ii)*abs(hr(n,j,ii))*sdev
                  if(sdev/abs(hr(n,j,ii)).lt.1.)then
		    dph=asin(sdev/abs(hr(n,j,ii)))/rad
                    write(18,999)deltat/freq(ii),rho,ph,rho+drho,
     $                           rho-drho,ph+dph,ph-dph,
     $                           fracthresh2(n,ii),fracdat2(n,ii)
                  else
                    write(18,999)deltat/freq(ii),rho,ph,rho+drho,
     $                           rho-drho,ph+180.,ph-180.,
     $                           fracthresh2(n,ii),fracdat2(n,ii)
                  endif
	        else
	          write(18,999)-999.,-999.,-999.,-999.,-999.,
     $                         -999.,-999.,-999.,-999.
	        endif
  987           continue
        close(unit=18)
      endif
  996 format(3(a,i8))
  997 format(a,5i8)
  998 format(3(a,g15.7))
  999 format(1x,12g15.7)
 1000 format(1x,8g15.7,1x,i1)
 1001 format(1x,9g15.7,1x,i1)
      end
