      program reads10
c
      implicit real(a-h,k-z)
      integer imm,is1,is2,itime,imin,imax
      parameter (im=36000,icon=1) 
      character fname*16,c7*7,c1*1,fname2*15,ctime*2

      real var(8,im),var2(8,im),var3(8,im),ave(8),cov(4)
c
      if(icon.eq.1) goto 1111
      open(11,file='filename.txt',form='formatted',status='old')
c
      imin=36000
      do 1000 ifile=1,240
      read(11,'(a16)') fname
      fname2='S0_'//fname(5:16)
      open(12,file='Orig/'//fname,form='formatted',status='old')
      open(22,file=fname2,form='formatted',status='unknown')
      iskip=mod(ifile-1,24)*36000
      
      itimeb=-1
      do it=1,im
      read(12,'(a7,8(a1,f7.3))',end=9912) 
     * c7,(c1,var(i,it),i=1,8)
      call c2int(c7,imm,is1,is2)
      itime=imm*600+is1*10+is2
      if(itime.ne.itimeb+1)  
     *write(667,*) fname(5:16),it,itime,itimeb
      write(22,'(i6,x,a6,8(x,f7.3))') 
     * iskip+itime,c7(1:2)//c7(4:7),(var(i,it),i=1,8)      
      imax=it
      itimeb=itime
      enddo
 9912 close(12)
      if(imax.lt.imin) imin=imax
      close(22)
      write(6,*) ' file = ',fname,'  ',fname2,imax

 1000 continue   
      close(11)
      write(6,*) ' data min = ',imin
      if(icon.eq.0.and.imin.lt.36000) stop
 1111 continue
c-----------------------------------------
      open(11,file='filename2.txt',form='formatted',status='old')
      open(21,file='s10all.gad',form='unformatted'
     *       ,status='unknown',recl=im,access='direct')
      open(24,file='s10movave.gad',form='unformatted'
     *       ,status='unknown',recl=im,access='direct')
      open(23,file='s10cov.dat',form='formatted',status='unknown')
      irec=0

      do 2000 ifile=1,240
      if(mod(ifile,24).eq.1)
     *write(6,'(3a)') '    file name    wT      wq     wT2     wq2  '
     *,'   u       v       w       T    '
     *,'   c       q       T2      q2    roua'   

c
      read(11,'(a15)') fname2
c      if(ifile.lt.130.or.ifile.gt.136) then
c      irec=irec+8
c      goto 2000
c      endif 
      open(12,file=fname2,form='formatted',status='old')
      do it=1,im
      read(12,'(i6,x,a6,8(x,f7.3))') 
     * itime,xtime,(var(i,it),i=1,8)      
      enddo
 9913 close(12)
c
c++++ removing spike noise ++++
      do i=1,4       
      aaa=0.2
      do it=2,im-1
      if(abs(var(i,it)-var(i,it-1)).gt.aaa.and.
     *   abs(var(i,it)-var(i,it+1)).gt.aaa) 
     * var(i,it)=(var(i,it-1)+var(i,it+1))*0.5  
      enddo
      enddo
c
      call rh2q(var,im,8,em,tm)
      call avevar(var,var2,ave,im,8)
      call movave(var2,var3,im,8,300)
c
c++++ density of air (kg/m3) ++++
      roua=(((1000.-em)/2.8705)+(em/4.61495))/tm

      call covar(var2,cov,im,8,roua)
c
      write(23,'(i3,x,a8,x,13g12.4)') 
     * ifile,fname2(4:11),(cov(i),i=1,4),(ave(i),i=1,8),roua
      write(6,'(i3,x,a8,x,13f8.2)') 
     * ifile,fname2(4:11),(cov(i),i=1,4),(ave(i),i=1,8),roua
c
      do i=1,8
      irec=irec+1
      write(21,rec=irec) (var2(i,it),it=1,im)
      write(24,rec=irec) (var3(i,it),it=1,im)
      enddo
 2000 continue
      close(11)
      close(21)
      close(23)
      close(24)

      stop
      end
c=============================================
      subroutine rh2q(var,im,jm,em,tm)
c------------------------------------
c var    u  v  w  T  H2O CO2 Tair RH
c------------------------------------
      integer im,jm,i
      real var(jm,im)
c++++ definition of function (trx: degC)++++
       fnce(trx)=6.1078d0*(10.d0**(7.5d0*trx/(237.3d0+trx)))
c

      em=0.
      tm=0.
      do i=1,im
      eee=fnce(var(7,i))*var(8,i)*0.01
      var(5,i)=40.*var(5,i) ! H2O (g/kg)
      var(6,i)=2.*var(6,i) ! CO2 (ppm)
      var(8,i)=1000.*0.622*eee/(1000.-0.378*eee) ! q (g/kg)
      em=em+eee
      tm=tm+var(7,i)
      enddo
      em=em/real(im)
      tm=273.15+tm/real(im)
c
      return
      end   
c=============================================
      subroutine avevar(var,var2,ave,im,jm)
c
      integer im,jm
      real var(jm,im),var2(jm,im),ave(jm)
c
      do j=1,jm
      sum=0.
      do i=1,im 
      sum=sum+var(j,i)
      enddo
      ave(j)=sum/real(im)
      do i=1,im
c      if(j.le.4.or.j.ge.7) then
      var2(j,i)=var(j,i)-ave(j)
c      else
c      var2(j,i)=var(j,i)
c      endif
      enddo
      enddo
c
      return
      end
c==============================================
      subroutine movave(var2,var3,im,jm,iave)
c
      integer im,jm,iave
      real var2(jm,im),var3(jm,im)
c
      do j=1,jm
      do i=1,im 
      sum=0.
      do ii=i-iave,i+iave
      iii=ii
      if(iii.lt.1) iii=ii+im
      if(iii.gt.im) iii=ii-im
      sum=sum+var2(j,iii)
      enddo
      var3(j,i)=sum/real(2*iave+1)
      if(j.eq.8) then
      ave=0.
      do ii=i-6000,i
      iii=ii
      if(iii.lt.1) iii=ii+im
      ave=ave+var2(j,iii)
      enddo
      ave=ave/6001.
      var3(j,i)=var3(j,i)-ave
      endif

      enddo
      enddo
c 
      return
      end
c================================================
      subroutine covar(var2,cov,im,jm,roua)
c
      integer im,jm
      real var2(jm,im),cov(4)
c--------------------------------------
c var2   u' v' w' T' c' q' T2' q2'
c cov    w'T' w'q' w'T2' w'q2' 
c--------------------------------------
c++++ specific heat of air (J/kgK) ++++
          cp = 1.01e3
c++++ latent heat of vaporization (J/kg)
      lamuda  = 2.50e6
c
      sum1=0. 
      sum2=0. 
      sum3=0. 
      sum4=0. 
      do i=1,im
      sum1=sum1+var2(3,i)*var2(4,i)
      sum2=sum2+var2(3,i)*var2(5,i)
      sum3=sum3+var2(3,i)*var2(7,i)
      sum4=sum4+var2(3,i)*var2(8,i)
      enddo
      cov(1)=roua*cp*sum1/real(im)
      cov(2)=1.e-3*roua*lamuda*sum2/real(im)
      cov(3)=roua*cp*sum3/real(im)
      cov(4)=1.e-3*roua*lamuda*sum4/real(im)
c
      return
      end
c=====================================
      subroutine c2int(c7,imm,is1,is2)
c
      integer imm,is1,is2,i1,i2,i3,i4,i5
      character c7*7,c1*1,c2*1,c3*1,c4*1,c5*1
c
      c1=c7(1:1)     
      c2=c7(2:2)     
      c3=c7(4:4)     
      c4=c7(5:5)     
      c5=c7(7:7)
c
      call c2int2(c1,i1)     
      call c2int2(c2,i2)     
      call c2int2(c3,i3)     
      call c2int2(c4,i4)     
      call c2int2(c5,i5)
c
      imm=i1*10+i2
      is1=i3*10+i4
      is2=i5           
c
      return
      end
c=======================
      subroutine c2int2(c1,i1)
      character c1*1
      integer i1
c
      if(c1.eq.'1') i1=1
      if(c1.eq.'2') i1=2
      if(c1.eq.'3') i1=3
      if(c1.eq.'4') i1=4
      if(c1.eq.'5') i1=5
      if(c1.eq.'6') i1=6
      if(c1.eq.'7') i1=7
      if(c1.eq.'8') i1=8
      if(c1.eq.'9') i1=9
      if(c1.eq.'0') i1=0
c
      return
      end
c===============================



