      subroutine bcoher2_b(nsect,ninp,nref,c2thresh,
     $                     ffts,c2flag,c2,info)
c
c birrp
c
c multiple coherence estimation on each section by block averaging
c
c  input variables:
c     nsect--number of sections in ffts
c     ninp--number of output (i.e., local magnetic field) time series
c     nref--number of input (i.e., remote magnetic field) time series
c     c2thresh-coherence threshold
c     ffts--complex array containing nsect section fft's
c           for ninp+nref variables
c  output variables:
c     c2flag-logical vector, .true. if coherence threshold is exceeded
c     c2--array containing nsect/nblk section average coherence
c         estimates for ninp output variables
c     info--error flag, 0 for normal return, 1 for singular matrix
c
      include 'parameters.h'
      parameter (nblk=10)
      complex ffts
      logical c2flag
      dimension ffts(nsectm,*),
     $          c2(nsectm/nblk,*),c2flag(*)
      double precision syy,rcond
      double complex sxx,sxy,syx,z,sypyp
      dimension sxx(nrefm,nrefm),syx(nrefm),sxy(nrefm),
     $          ipvt(nrefm),z(nrefm)
      info=0
      do 50 n=0,nsect-1,nblk
        do 40 m=1,ninp
          do 10 i=1,nref
            syx(i)=0.
            do 10 j=1,nref
   10         sxx(i,j)=0.
          syy=0.
          do 20 ii=n+1,n+nblk
            syy=syy+abs(ffts(ii,m))**2
            do 20 i=1,nref
              syx(i)=syx(i)+conjg(ffts(ii,m))*ffts(ii,i+ninp)
              do 20 j=1,nref
   20           sxx(i,j)=sxx(i,j)+conjg(ffts(ii,i+ninp))*
     $                   ffts(ii,j+ninp)
          do 25 i=1,nref
   25       sxy(i)=conjg(syx(i))
          call zgeco(sxx,nrefm,nref,ipvt,rcond,z)
          if(1.d0+rcond.eq.1.d0)info=1
          job=0
          call zgesl(sxx,nrefm,nref,ipvt,sxy,job)
          sypyp=0. 
          do 30 i=1,nref
   30       sypyp=sypyp+syx(i)*sxy(i)
   40     c2(n/nblk+1,m)=dreal(sypyp)/syy
   50   continue
      if(c2thresh.eq.0.)return
      do 80 nn=1,ninp
        do 80 i=1,nsect/nblk
          if(c2(i,nn).le.c2thresh)c2flag(i)=.false.
   80     continue
      return
      end
      subroutine bcoher2_e(nsect,n1,nout,ninp,c2thresh,
     $                     ffts,c2flag,c2,info)
c
c birrp
c
c multiple coherence estimation on each section by block averaging for remote
c reference data
c
c  input variables:
c     nsect--number of sections in ffts
c     n1--number of output time series to process
c     nout--number of output (i.e., electric field) time series
c     ninp--number of input (i.e., magnetic field)time series
c     c2thresh-coherence threshold
c     ffts--complex array containing nsect section fft's
c           for nout+ninp variables
c  output variables:
c     c2flag-logical vector, .true. if coherence threshold is exceeded
c     c2--array containing nsect/nblk section average coherence
c         estimates for nout output variables
c     info--error flag, 0 for normal return, 1 for singular matrix
c
      include 'parameters.h'
      parameter (nblk=10)
      complex ffts
      logical c2flag
      dimension ffts(nsectm,*),
     $          c2(nsectm/nblk,*),c2flag(*)
      double precision syy,rcond
      double complex sxrx,sxry,sxx,syx,a,sypyp,syyp
      dimension sxrx(ninpm,ninpm),sxry(ninpm),
     $          sxx(ninpm,ninpm),syx(ninpm),
     $          ipvt(ninpm),a(ninpm)
      info=0
      do 50 n=0,nsect-1,nblk
        do 40 m=1,n1
          do 10 i=1,ninp
            sxry(i)=0.
            syx(i)=0.
            do 10 j=1,ninp
              sxx(i,j)=0.
   10         sxrx(i,j)=0.
          syy=0.
          do 20 ii=n+1,n+nblk
            syy=syy+abs(ffts(ii,m))**2
            do 20 i=1,ninp
              sxry(i)=sxry(i)+conjg(ffts(ii,i+nout+ninp))*
     $             ffts(ii,m)
              syx(i)=syx(i)+conjg(ffts(ii,m))*ffts(ii,i+nout)
              do 20 j=1,ninp
                sxx(i,j)=sxx(i,j)+conjg(ffts(ii,i+nout))*
     $                   ffts(ii,j+nout)
   20           sxrx(i,j)=sxrx(i,j)+conjg(ffts(ii,i+nout+ninp))*
     $                 ffts(ii,j+nout)
          call zgeco(sxrx,ninpm,ninp,ipvt,rcond,a)
          if(1.d0+rcond.eq.1.d0)info=1
          job=0
          call zgesl(sxrx,ninpm,ninp,ipvt,sxry,job)
          sypyp=0.
          syyp=0.
          do 30 i=1,ninp
            a(i)=0.
            do 25 j=1,ninp
   25         a(i)=a(i)+sxx(i,j)*sxry(j)
            syyp=syyp+syx(i)*sxry(i)
   30       sypyp=sypyp+conjg(sxry(i))*a(i)
   40     c2(n/nblk+1,m)=abs(syyp)**2/(syy*dreal(sypyp))
   50   continue
      if(c2thresh.eq.0.)return
      do 80 i=1,nsect/nblk
        do 80 nn=1,n1
          if(c2(i,nn).le.c2thresh)c2flag(i)=.false.
   80     continue
      return
      end
      subroutine ccoher1(nsct,ninp,nref,nn,ftsec,h,hd1,wt,gd1,g,
     $                  glo,ghi)
c
c  birrp
c
c  compute complex coherence between measured and predicted output variables
c
c  input variables:
c  
c  nsct--number of sections
c  ninp--number of output variables
c  nref--number of input variables
c  nn--index of output variable
c  ftsec--complex array containing nsct ffts
c  h--complex array containing transfer functions between ninp output and
c      nref input variables
c  hd1--complex array of delete one transfer function estimates
c  wt--weight vector
c  gd1--work array of length at least nsct
c  
c  output variables
c  g--complex coherence (not squared coherence!)
c  glo--lower 95% confidence limit on coherence magnitude
c  ghi--upper 95% confidence limit on coherence magnitude
      include 'parameters.h'
      complex ftsec,h,hd1,g
      dimension ftsec(nsectm,*),h(ninpm,*),
     $          hd1(nrefm+1,*),wt(*),gd1(*)
      double precision pwr,gd1av,syy,var
      double complex bb,sxx,yps,zt0,t1,sxx1,yp
      dimension bb(nrefm),sxx(nrefm,nrefm),yp(nrefm)
      pwr=0.
      do 10 i=1,nref
        bb(i)=0.
        do 10 j=1,nref
   10     sxx(i,j)=0.
      do 20 ii=1,nsct
        pwr=pwr+(wt(ii)*abs(ftsec(ii,nn)))**2
        do 20 i=1,nref
           bb(i)=bb(i)+wt(ii)*wt(ii)*ftsec(ii,ninp+i)*
     $                 conjg(ftsec(ii,nn))
           do 20 j=1,nref
   20        sxx(i,j)=sxx(i,j)+wt(ii)*wt(ii)*
     $                conjg(ftsec(ii,ninp+i))*
     $                ftsec(ii,ninp+j)
      yps=0.
      zt0=0.
      do 30 i=1,nref
        yp(i)=0.
        do 25 j=1,nref
   25     yp(i)=yp(i)+sxx(i,j)*h(nn,j)
        yps=yps+conjg(h(nn,i))*yp(i)
   30   zt0=zt0+bb(i)*h(nn,i)
      g=zt0/sqrt(pwr*yps)
c  compute jackknifed estimate of coherence 95% confidence interval
      gd1av=0.
      nd=0
      do 50 ii=1,nsct
        if(wt(ii).gt.0.)then
          nd=nd+1
          syy=pwr-(wt(ii)*abs(ftsec(ii,nn)))**2
          zt0=0.
          yps=0.
          do 40 i=1,nref
            t1=bb(i)-wt(ii)*wt(ii)*ftsec(ii,ninp+i)*
     $               conjg(ftsec(ii,nn))
            zt0=zt0+t1*hd1(i,ii)
            yp(i)=0.
            do 35 j=1,nref
              sxx1=sxx(i,j)-wt(ii)*wt(ii)*conjg(ftsec(ii,ninp+i))*
     $             ftsec(ii,ninp+j)
   35         yp(i)=yp(i)+sxx1*hd1(j,ii)
   40       yps=yps+conjg(hd1(i,ii))*yp(i)
            gd=abs(zt0/sqrt(syy*yps))
            gd1(nd)=arctanh(gd)
            gd1av=gd1av+gd1(nd)
        endif
   50   continue
      gd1av=gd1av/nd
      var=0.
      do 60 i=1,nd
   60   var=var+(gd1av-gd1(i))**2
      err=sqrt((nd-1)*var/nd)
      g1=arctanh(abs(g))
      glo=max(0.,tanh1(g1-1.96*err))
      ghi=min(tanh1(g1+1.96*err),1.)
      return
      end
      subroutine ccoher2(nsct,nout,ninp,nn,ftsec,h,hd1,wt,gd1,g,
     $                  glo,ghi)
c
c  birrp
c
c  compute complex coherence between measured and predicted output variables
c
c  input variables:
c  
c  nsct--number of sections
c  nout--number of output variables
c  ninp--number of input variables
c  nn--index of output variable
c  ftsec--complex array containing nsct ffts
c  h--complex array containing transfer functions between nout output and
c      ninp input variables
c  hd1--complex array of delete one transfer function estimates
c  wt--weight vector
c  gd1--work array of length at least nsct
c  
c  output variables
c  g--complex coherence (not squared coherence!)
c  glo--lower 95% confidence limit on coherence magnitude
c  ghi--upper 95% confidence limit on coherence magnitude
      include 'parameters.h'
      complex ftsec,h,hd1,g
      dimension ftsec(nsectm,*),h(noutm,*),
     $          hd1(nrefm+1,*),wt(*),gd1(*)
      double precision pwr,gd1av,syy,var
      double complex bb,sxx,yps,zt0,t1,sxx1,yp
      dimension bb(ninpm),sxx(ninpm,ninpm),yp(ninpm)
      pwr=0.
      do 10 i=1,ninp
        bb(i)=0.
        do 10 j=1,ninp
   10     sxx(i,j)=0.
      do 20 ii=1,nsct
        pwr=pwr+(wt(ii)*abs(ftsec(ii,nn)))**2
        do 20 i=1,ninp
           bb(i)=bb(i)+wt(ii)*wt(ii)*ftsec(ii,nout+i)*
     $                 conjg(ftsec(ii,nn))
           do 20 j=1,ninp
   20        sxx(i,j)=sxx(i,j)+wt(ii)*wt(ii)*
     $                conjg(ftsec(ii,nout+i))*
     $                ftsec(ii,nout+j)
      yps=0.
      zt0=0.
      do 30 i=1,ninp
        yp(i)=0.
        do 25 j=1,ninp
   25     yp(i)=yp(i)+sxx(i,j)*h(nn,j)
        yps=yps+conjg(h(nn,i))*yp(i)
   30   zt0=zt0+bb(i)*h(nn,i)
      g=zt0/sqrt(pwr*yps)
c  compute jackknifed estimate of coherence 95% confidence interval
      gd1av=0.
      nd=0
      do 50 ii=1,nsct
        if(wt(ii).gt.0.)then
          nd=nd+1
          syy=pwr-(wt(ii)*abs(ftsec(ii,nn)))**2
          zt0=0.
          yps=0.
          do 40 i=1,ninp
            t1=bb(i)-wt(ii)*wt(ii)*ftsec(ii,nout+i)*
     $               conjg(ftsec(ii,nn))
            zt0=zt0+t1*hd1(i,ii)
            yp(i)=0.
            do 35 j=1,ninp
              sxx1=sxx(i,j)-wt(ii)*wt(ii)*conjg(ftsec(ii,nout+i))*
     $             ftsec(ii,nout+j)
   35         yp(i)=yp(i)+sxx1*hd1(j,ii)
   40       yps=yps+conjg(hd1(i,ii))*yp(i)
            gd=abs(zt0/sqrt(syy*yps))
            gd1(nd)=arctanh(gd)
            gd1av=gd1av+gd1(nd)
        endif
   50   continue
      gd1av=gd1av/nd
      var=0.
      do 60 i=1,nd
   60   var=var+(gd1av-gd1(i))**2
      err=sqrt((nd-1)*var/nd)
      g1=arctanh(abs(g))
      if(g1-1.96*err.lt.0.)then
        glo=0.
      else
        glo=tanh1(g1-1.96*err)
      endif
      ghi=min(tanh1(g1+1.96*err),1.)
      return
      end
      real function arctanh(z)
      arctanh=0.5*log((1.+z)/(1.-z))
      return
      end
      real function tanh1(z)
      tanh1=(1.-exp(-2.*z))/(1.+exp(-2.*z))
      return
      end

