      program wind_histo
      implicit none

      integer*4 imxrec,imxbu,imxbd
      parameter(imxrec=1000000,imxbu=20,imxbd=360)
      integer*4 i,nrec,nbin_u,nbin_d,j,k,nu(imxbu,imxbd),nd(imxbd)
      integer*4 isbad(imxrec),nrec1,neq
      real*8 r1,r2,r3,r4,r5,r6,r7,r8,u(imxrec),d(imxrec),ds,u1(imxrec),
     &       d1(imxrec),du
      real*8 bd(imxbd),bu(imxbu),avg_u(imxbu)
      character finp*200

! reads the 1st argument of the calling command, the
! name of the input ASCII meteo archive:
      call getarg(1,finp) 
      open(1,file=finp,status='old')
        read(1,*) ! skip file header
        do i=1,imxrec
          read(1,*,end=108) r1,r2,r3,r4,r5,r6,r7,r8,u1(i),d1(i)
c          write(*,*) u1(i),d1(i)
        end do
        stop'FERR: imxrec too small !!'
  108 close(1)
      nrec1 = i - 1 ! N. of records in the archive
      write(*,*)'N. of records in archive = ',nrec1

! ================ Filtering of data =================
! First we initialize to "good" (isbad(i)=0 all archive records
      do i=1,imxrec
        isbad(i) = 0
      end do

! 1) Checks whether data are in physically meaningful ranges
! Must be: 0 < = u1(i) < = 25 m/s AND 0 < = d1(i) < = 360
      do i=1, nrec1
        if(u1(i).lt.0.d0 .OR. u1(i).gt. 25.d0) isbad(i)=1 
! .... d1(i) ?
      end do

! 2) Wind speed should change more than du m/s in 12 hours
! if records are every 30 min ==> 12 hrs = 24 records
! du = wind speed resolution in m/s
      du = 0.1d0 ! to be modified ...
      do i=24,nrec1
        neq = 1
        do j = i, i-24+1, -1
         if(dabs(u1(i)-u1(j)).ge.du) neq = 0
        end do
        if (neq.ne.0) then
          do j = i, i-24+1, -1
            isbad(j) = 1
          end do
        end if
      end do
! 3) Additional filters to be added ???
! ......

! Final filter step: copy good data into new vectors u(i) and d(i):
      nrec = 0
      do i=1,nrec1
        if (isbad(i).eq.0) then
          nrec = nrec + 1
          u(nrec) = u1(i)
          d(nrec) = d1(i)
        end if
      end do
      write(*,*)'N. of records after filtering: ',nrec
! ================ End filtering of data =================


! === Definition of direction bins (the edges):
      nbin_d = 17
      ds = 360.d0 / dble(nbin_d-1)
      do j=1,nbin_d+1
        bd(j) = ds*(j-1) - ds*.5d0
        nd(j) = 0
      end do
      bd(1) = 0.d0
      bd(nbin_d+1) = 360.d0

! write out bins, uncomment to debug:     
!      do j=1,nbin_d+1
!        write(*,*) j, bd(j)
!      end do

! Definition of the velocity bins (the edges):
      nbin_u = 6
      bu(1)=0.d0; bu(2)=1.d0; bu(3)=2.d0 
      bu(4)=3.d0; bu(5)=4.d0; bu(6)=5.d0; bu(7)=100.d0
      do k=1,nbin_u
        do j=1,nbin_d        
          nu(k,j) = 0
        end do
      end do

! calculates the population of the bins:
      do j=1,nbin_d
        avg_u(j)=0.d0
        do i=1,nrec
          if (bd(j).le.d(i).AND.d(i).lt.bd(j+1))then
            nd(j) = nd(j) + 1
            avg_u(j) = avg_u(j) + u(i)
            do k=1,nbin_u
              if(bu(k).le.u(i).AND.u(i).lt.bu(k+1)) nu(k,j) = nu(k,j)+1
            end do
          end if
        end do ! end loop on archive records "i"
      end do ! end loop on direction bins "j"
 
! the first and the last direction bins (around the North) are merged: 
      nd(1) = nd(1)+nd(nbin_d)
      do k=1,nbin_u
        nu(k,1) = nu(k,1) + nu(k,nbin_d)
      end do

! calculation of the average velocity in each direction bin:
      do j=1,nbin_d - 1
        avg_u(j) = avg_u(j) / dble(nd(j))
      end do


! writing outputs:
      open(2,file='frequencies.dat',status='unknown')
      open(3,file='bin_populat.dat',status='unknown')
      write(2,'(a)')
     &  '#  DIREZ   Vdire   V_0-1   V_1-2   V_2-3   '//
     &  'V_3-4   V_4-5    V>=5'
      write(3,'(a)')
     &  '#  DIREZ   Ndire   N_0-1   N_1-2   N_2-3   '//
     &  'N_3-4   N_4-5    N>=5   AVG_U(m/s)'
      ds = 360.d0 / (nbin_d-1)
      do j=1,nbin_d-1
        r1 = ds * (j-1)
        write(2,10)r1,nd(j)/dble(max(nrec,1)),
     &       (nu(k,j)/dble(max(nd(j),1)),k=1,nbin_u)
        write(3,20)r1,nd(j),(nu(k,j),k=1,nbin_u),avg_u(j)
      end do
      close(2)
      close(3)

10    FORMAT(F8.1,10F8.3)
20    FORMAT(F8.1,7I8,3F8.3)

      stop'INFO: ==== NORMAL COMPLETION ===='
      end
