* fft_2hr_all(6ch data)  with linear trand erased 
*
      integer tim(5),tim0(5),nmax,np,nf,nf2,nb,nb2,nb4
      real ux(72000),uy(72000),uz(72000)
      real u(72000),v(72000),w(72000)
      real t(72000),q(72000),c(72000)
      real dat(6),adat(6),dt,a1,a2,au,av,aw
      real s(16384),frq(32768),dbt
      real wind_dir,hpass
      real ww,uu,vv,tt,qq,cc
      real wu,wv,wt,wq,wc
* fourie parameter
      complex fw(32768),fu(32768),fv(32768)
      complex ft(32768),fq(32768),fc(32768)
      real fwr(32768),fur(32768),fvr(32768)
      real ftr(32768),fqr(32768),fcr(32768)
      real fwi(32768),fui(32768),fvi(32768)
      real fti(32768),fqi(32768),fci(32768)
* power spectra (single)
      real pw(32768),pu(32768),pv(32768)
      real pt(32768),pq(32768),pc(32768)
      real pw1(32768),pu1(32768),pv1(32768)
      real pt1(32768),pq1(32768),pc1(32768)
* co-spectra
      complex pwu(32768),pwv(32768)
      complex pwt(32768),pwq(32768),pwc(32768)
      complex pwu1(32768),pwv1(32768)
      complex pwt1(32768),pwq1(32768),pwc1(32768)
      real pwur(32768),pwvr(32768)
      real pwtr(32768),pwqr(32768),pwcr(32768)
      real pwur1(32768),pwvr1(32768)
      real pwtr1(32768),pwqr1(32768),pwcr1(32768)
      real pwui(32768),pwvi(32768)
      real pwti(32768),pwqi(32768),pwci(32768)
      real pwui1(32768),pwvi1(32768)
      real pwti1(32768),pwqi1(32768),pwci1(32768)
* paramters
      dt=0.1
      nf=50
      nf2=nf/2
*
* data reading
      do 100 i=0,999999,1
        nmax=0
        do 101 k=1,6,1
          adat(k)=0.
101     continue
        do 200 j=1,72000,1
          read(*,*,end=999) tim,dat
          if(j.eq.1) then
            do 201 k=1,5,1
              tim0(k)=tim(k)
201         continue
          end if
          do 202 k=1,3,1
            if(dat(k).lt.-99.) goto 200
202       continue
          dat(5)=dat(5)*44.0
          dat(6)=dat(6)*18.0/1000.
          do 203 k=1,6,1
            adat(k)=adat(k)+dat(k)
203       continue
          nmax=nmax+1
          ux(nmax)=dat(1)
          uy(nmax)=dat(2)
          uz(nmax)=dat(3)
          t(nmax)=dat(4)
          q(nmax)=dat(6)
          c(nmax)=dat(5)
200     continue
        if(nmax.lt.66000) goto 100
*
*        do 205 j=1,16500,1
*          jj=(j-1)*4
*          ux(j)=(ux(jj+1)+ux(jj+2)+ux(jj+3)+ux(jj+4))/4.
*          uy(j)=(uy(jj+1)+uy(jj+2)+uy(jj+3)+uy(jj+4))/4.
*          uz(j)=(uz(jj+1)+uz(jj+2)+uz(jj+3)+uz(jj+4))/4.
*          t(j)=(t(jj+1)+t(jj+2)+t(jj+3)+t(jj+4))/4.
*          q(j)=(q(jj+1)+q(jj+2)+q(jj+3)+q(jj+4))/4.
*          c(j)=(c(jj+1)+c(jj+2)+c(jj+3)+c(jj+4))/4.
*205     continue
*        nmax=nmax/4
*
        do 102 k=1,6,1
          adat(k)=adat(k)/real(nmax)
102     continue
        if(adat(1).gt.0.) then
          wind_dir=atan(-adat(2)/adat(1))*180./3.141592
        else
          wind_dir=180+atan(-adat(2)/adat(1))*180./3.141592
        end if
        wind_dir=wind_dir-30+10+180
        if(wind_dir.lt.0.) wind_dir=wind_dir+360.
        if(wind_dir.ge.360.) wind_dir=wind_dir-360.
        a1=adat(1)*adat(1)+adat(2)*adat(2)
        a2=a1+adat(3)*adat(3)
        if(a1.gt.0.) a1=a1**0.5
        if(a2.gt.0.) a2=a2**0.5
        au=0.
        av=0.
        aw=0.
        do 300 j=1,nmax,1
          u(j)=(adat(1)*ux(j)+adat(2)*uy(j))/a1
          v(j)=(adat(1)*uy(j)-adat(2)*ux(j))/a1
          w(j)=(a1*uz(j)-adat(3)*u(j))/a2
          au=au+u(j)/real(nmax)
          av=av+v(j)/real(nmax)
          aw=aw+w(j)/real(nmax)
300     continue
        do 301 j=1,nmax,1
          ux(j)=ux(j)-adat(1)
          uy(j)=uy(j)-adat(2)
          uz(j)=uz(j)-adat(3)
          t(j)=t(j)-adat(4)
          q(j)=q(j)-adat(6)
          c(j)=c(j)-adat(5)
          u(j)=u(j)-au
          v(j)=v(j)-av
          w(j)=w(j)-aw
301     continue
*
* high pass
*	hpass=2.
*        do 310 j=1,nmax-int(hpass)+1,1
*          do 311 k=1,int(hpass)-1,1
*	    u(j)=u(j)+u(j+k)
*	    v(j)=v(j)+v(j+k)
*	    w(j)=w(j)+w(j+k)
*            t(j)=t(j)+t(j+k)
*            q(j)=q(j)+q(j+k)
*            c(j)=c(j)+c(j+k)
*            ux(j)=ux(j)+ux(j+k)
*            uy(j)=uy(j)+uy(j+k)
*            uz(j)=uz(j)+uz(j+k)
*311       continue
*          u(j)=u(j)/hpass
*          v(j)=v(j)/hpass
*          w(j)=w(j)/hpass
*          t(j)=t(j)/hpass
*          q(j)=q(j)/hpass
*          c(j)=c(j)/hpass
*          ux(j)=ux(j)/hpass
*          uy(j)=uy(j)/hpass
*          uz(j)=uz(j)/hpass
*310     continue
*
* remove linear trend of T, Q, C
        var_t=0.
        var_q=0.
        var_c=0.
        cov_t=0.
        cov_q=0.
        cov_c=0.
        avg_j=real(nmax-1)/2.
        do 350 j=1,nmax,1
          var_t=var_t+t(j)*t(j)/real(nmax)
          var_q=var_q+q(j)*q(j)/real(nmax)
          var_c=var_c+c(j)*c(j)/real(nmax)
          cov_t=cov_t+t(j)*real(j-1)/real(nmax)
          cov_q=cov_q+q(j)*real(j-1)/real(nmax)
          cov_c=cov_c+c(j)*real(j-1)/real(nmax)
350     continue
        do 351 j=1,nmax,1
          t(j)=t(j)-real(j-1)*cov_t/var_t
          q(j)=q(j)-real(j-1)*cov_q/var_q
          c(j)=c(j)-real(j-1)*cov_c/var_c
351     continue
*
* variance
        ww=0.
        uu=0.
        vv=0.
        tt=0.
        qq=0.
        cc=0.
        wu=0.
        wv=0.
        wt=0.
        wq=0.
        wc=0.
        do 600 j=1,nmax,1
          ww=ww+w(j)*w(j)/real(nmax)
          uu=uu+u(j)*u(j)/real(nmax)
          vv=vv+v(j)*v(j)/real(nmax)
          tt=tt+t(j)*t(j)/real(nmax)
          qq=qq+q(j)*q(j)/real(nmax)
          cc=cc+c(j)*c(j)/real(nmax)
          wu=wu+w(j)*u(j)/real(nmax)
          wv=wv+w(j)*v(j)/real(nmax)
          wt=wt+w(j)*t(j)/real(nmax)
          wq=wq+w(j)*q(j)/real(nmax)
          wc=wc+w(j)*c(j)/real(nmax)
600     continue
*
* fourie transform
        np=alog(real(nmax)+0.5)/alog(2.0)
        nb=2**np
        nb2=nb/2
        nb4=nb/4
        do 400 j=1,nb4-1,1
          s(j)=sin(real(j)*6.283185/real(nb))
400     continue
        call fft(w,fwr,fwi,nb,nb2,nb4,np,s)
        call fft(u,fur,fui,nb,nb2,nb4,np,s)
        call fft(v,fvr,fvi,nb,nb2,nb4,np,s)
        call fft(t,ftr,fti,nb,nb2,nb4,np,s)
        call fft(q,fqr,fqi,nb,nb2,nb4,np,s)
        call fft(c,fcr,fci,nb,nb2,nb4,np,s)
        fwi(1)=0.
        fui(1)=0.
        fvi(1)=0.
        fti(1)=0.
        fqi(1)=0.
        fci(1)=0.
        do 410 j=1,nb2,1
          fw(j)=fwr(j)+fwi(j)*(0.0,1.0)
          fu(j)=fur(j)+fui(j)*(0.0,1.0)
          fv(j)=fvr(j)+fvi(j)*(0.0,1.0)
          ft(j)=ftr(j)+fti(j)*(0.0,1.0)
          fq(j)=fqr(j)+fqi(j)*(0.0,1.0)
          fc(j)=fcr(j)+fci(j)*(0.0,1.0)
410     continue
*
* spectrum
        tn=2.0*dt/real(nb-1)
        dbt=1.0/real(nb-1)/dt*real(nf2)
        do 500 j=1,nb2,1
          pw1(j)=fw(j)*conjg(fw(j))*tn
          pu1(j)=fu(j)*conjg(fu(j))*tn
          pv1(j)=fv(j)*conjg(fv(j))*tn
          pt1(j)=ft(j)*conjg(ft(j))*tn
          pq1(j)=fq(j)*conjg(fq(j))*tn
          pc1(j)=fc(j)*conjg(fc(j))*tn
          pwu1(j)=fw(j)*conjg(fu(j))*tn
          pwv1(j)=fw(j)*conjg(fv(j))*tn
          pwt1(j)=fw(j)*conjg(ft(j))*tn
          pwq1(j)=fw(j)*conjg(fq(j))*tn
          pwc1(j)=fw(j)*conjg(fc(j))*tn
          pwur1(j)=real(pwu1(j))
          pwvr1(j)=real(pwv1(j))
          pwtr1(j)=real(pwt1(j))
          pwqr1(j)=real(pwq1(j))
          pwcr1(j)=real(pwc1(j))
          pwui1(j)=aimag(pwu1(j))
          pwvi1(j)=aimag(pwv1(j))
          pwti1(j)=aimag(pwt1(j))
          pwqi1(j)=aimag(pwq1(j))
          pwci1(j)=aimag(pwc1(j))
500     continue
        call filter(pw1,pw,nb2,mf,nf,nf2)
        call filter(pu1,pu,nb2,mf,nf,nf2)
        call filter(pv1,pv,nb2,mf,nf,nf2)
        call filter(pt1,pt,nb2,mf,nf,nf2)
        call filter(pq1,pq,nb2,mf,nf,nf2)
        call filter(pc1,pc,nb2,mf,nf,nf2)
        call filter(pwur1,pwur,nb2,mf,nf,nf2)
        call filter(pwvr1,pwvr,nb2,mf,nf,nf2)
        call filter(pwtr1,pwtr,nb2,mf,nf,nf2)
        call filter(pwqr1,pwqr,nb2,mf,nf,nf2)
        call filter(pwcr1,pwcr,nb2,mf,nf,nf2)
        call filter(pwui1,pwui,nb2,mf,nf,nf2)
        call filter(pwvi1,pwvi,nb2,mf,nf,nf2)
        call filter(pwti1,pwti,nb2,mf,nf,nf2)
        call filter(pwqi1,pwqi,nb2,mf,nf,nf2)
        call filter(pwci1,pwci,nb2,mf,nf,nf2)
        do 510 j=1,nb,1
          frq(j)=dbt*real(j-1)
          pwu(j)=pwur(j)+pwui(j)*(0.0,1.0)
          pwv(j)=pwvr(j)+pwvi(j)*(0.0,1.0)
          pwt(j)=pwtr(j)+pwti(j)*(0.0,1.0)
          pwq(j)=pwqr(j)+pwqi(j)*(0.0,1.0)
          pwc(j)=pwcr(j)+pwci(j)*(0.0,1.0)
510     continue
*
* output
        open(31,file="avg.dat")
        open(30,file="spectr.dat")
        write(30,1031) tim0(2),tim0(3),tim0(4)
        write(30,1032) wind_dir,aw,au,av,adat(4),adat(6),adat(5)
        write(30,1033) ww,uu,vv,tt,qq,cc,wu,wv,wt,wq,wc
        write(31,1035) tim0(2),tim0(3),tim0(4),wind_dir,aw,au,av
     &                ,adat(4),adat(6),adat(5)
     &                ,ww,uu,vv,tt,qq,cc,wu,wv,wt,wq,wc 
	write(30,1000) 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.
        do 610 j=2,mf,1 
          write(30,1000) frq(j),frq(j)/au
     &          ,pw(j)/ww*frq(j)*dt
     &          ,pu(j)/uu*frq(j)*dt
     &          ,pv(j)/vv*frq(j)*dt
     &          ,pt(j)/tt*frq(j)*dt
     &          ,pq(j)/qq*frq(j)*dt
     &          ,pc(j)/cc*frq(j)*dt
     &          ,pwur(j)/wu*frq(j)*dt
     &          ,pwvr(j)/wv*frq(j)*dt
     &          ,pwtr(j)/wt*frq(j)*dt
     &          ,pwqr(j)/wq*frq(j)*dt
     &          ,pwcr(j)/wc*frq(j)*dt
610     continue
        close(30)
1000    format(f10.6,12(',',f10.6))
1031    format(3i6)
1032    format(10x,7(',',f10.6))
1033    format(10x,',',10x,11(',',f10.6))
1035    format(3i6,7(',',f10.6),11(',',f10.6))
*
        call mkname()
*
100   continue
999   continue
      end
*********************************************
      subroutine filter(ssp,sp,m,mf,nf,nf2)
      real ssp(32768),sp(32768)
      real sum
      integer m,mf,nf,nf2
      sum=0.
      mf=0
      do 100 k=1,m,nf2
	mf=mf+1
	sum=ssp(k)*nf
	if(nf.eq.1) goto 199
	do 150 i=1,nf-1,1
	  k1=k-i
	  k2=k+i
	  if(k1.le.0) k1=-(k1-1)+1
	  if(k2.ge.m) k2=m-(k2-m)
	  sum=sum+(ssp(k1)+ssp(k2))*real(nf-1)
150	continue
199	sp(mf)=sum/real(nf*nf)
100   continue
      end
*
      subroutine fft(x,a,b,n,n2,n4,np,s)
      real x(72000),a(32768),b(32768),s(16384)
      integer n,n2,n4,np
      integer ip,m1,m0,m2
      real s1,s2
      do 100 i=1,n2,1
	a(i)=x(i)+x(i+n2)
	b(i)=x(i)-x(i+n2)
100   continue
      do 200 i=1,n2,1
	x(2*i-1)=a(i)
	x(2*i)=b(i)
200   continue
      ip=n4
      m1=1
      do 300 m=1,np-1
	m2=2*m1
	m0=m1-1
	do 310 i=1,ip
	  i1=m2*(i-1)+1
	  i2=i1+m1
	  i3=i1+n2
	  i4=i3+m1
	  a(i1)=x(i1)+x(i3)
	  b(i1)=x(i1)-x(i3)
	  a(i2)=x(i2)
	  b(i2)=x(i4)
	  if(m.eq.1) goto 310
	  do 320 k=1,m0,1
	    s1=s(n4-ip*k)*x(i3+k)-s(ip*k)*x(i4+k)
	    s2=s(ip*k)*x(i3+k)+s(n4-ip*k)*x(i4+k)
	    i5=i1+k
	    i6=i1-k+m2
	    a(i5)=x(i5)+s1
	    a(i6)=x(i5)-s1
	    b(i5)=x(i2+k)+s2
	    b(i6)=-x(i2+k)+s2
320	  continue
310	continue
	do 350 j=1,ip,1
	  j1=m2*(j-1)
	  j2=2*j1
	  j3=j2+m2
	  do 360 k=1,m2
	    x(j2+k)=a(j1+k)
	    x(j3+k)=b(j1+k)
360	  continue
350	continue
	if(ip.eq.1) goto 999
	m1=m2
	ip=ip/2
300   continue
999   end
	
      
