!$Author$
!$Date$
!$Revision$
!$HeadURL$

!          INTERP_WDB:  weighted interpolation of wind data records
!                         03/20/09

!   Takes individual wdb files and weighting coefficients on the command
! line and does a weighted average of the wdb file values, putting the
! result into a new wdb file (also specified on the command line)

! Call by invoking:

! interp_wdb [new-file-name] [file-1 coeff-1] [file-2 coeff-2] ... [file-n coeff-n]

! Coefficients will be checked for summing to one.
! If not, a warning will be generated and the coefficieints normalized
! so that they sum to one.

      integer iargc
      character*50 argv
      integer idx, jdx, kdx, ndx, nfile, idxm1, idxs
      character*50 outfile, filename
      character*50, allocatable :: infile(:)
      real, allocatable :: weight(:)
      real wsum

      ! parameters for data records
      integer nbin, ibin
      parameter (nbin = 25) ! there are 27 bins, but only 25 recorded in a file
      real latitude, longitude
      real dirf(12,16), calm(12), maxminr(12), hrmax(12)
      real freq(12,16,nbin) !cumulative freq. < cutoff wind speed (cdf)
      integer lat_deg, lat_min, lon_deg, lon_min, elev
      character*40 pound, wmoid, name, country, state
      character lat_NS, lon_EW

      ! parameters for new weighted average data records
      real, allocatable :: trx(:), try(:), trz(:)
      real sc, px, py, pz, vnrm
      real wlatitude, wlongitude, welev
      real wcalm(12), wmaxminr(12), whrmax(12)
      integer iwhrmax(12)
      real wdirf(12,16), wfreq(12,16,nbin)

      ! parpmeter to convert degrees to radians
      sc = ATAN(1.)/45.

      ! get output file name and [file.wdb coefficient] pairs from command line
      if( iargc() .gt. 0 ) then
        if( mod(iargc()-1,2) .gt. 0 ) then
          write(*,*) 'file name coefficient pairs are required'
          stop
        end if
        nfile = (iargc()-1)/2
        allocate ( infile(1:nfile), weight(1:nfile) )
        allocate ( trx(1:nfile), try(1:nfile), trz(1:nfile) )
        do idx = 1, iargc()
          call getarg(idx,argv)
          if (idx .eq.1) then
            outfile = argv ! ouput file name
          else
            ! which part of pair
            idxm1 = idx - 1
            idxs = idx / 2
            if( mod(idxm1,2) .gt. 0 ) then
              ! input file name
              infile(idxs) = argv
            else
              ! weighting coefficient
              read(argv(1:),*) weight(idxs)
            end if
          end if
        end do
      else
        write(*,*) 'Arguments are required'
        write(*,*) 'interp_wdb [nw-f] [f-1 c-1] [f-2 c-2] ... [f-n c-n]'
        stop
      endif

      ! normalize weights to 1
      wsum = 0.0
      do ndx = 1, nfile
        wsum = wsum + weight(ndx)
      end do
      do ndx = 1, nfile
        weight(ndx) = weight(ndx) / wsum
      end do

      ! zero accumulator arrays
      welev = 0.0
      do idx = 1, 12
        wcalm(idx) = 0.0
        wmaxminr(idx) = 0.0
        whrmax(idx) = 0.0
        do jdx = 1, 16
          ! zero out direction arrays
          wdirf(idx,jdx) = 0.0
          do kdx = 1, nbin
            ! set histogram arrays so end is filled where not read
            wfreq(idx,jdx,kdx)=0.0
          end do
        end do
      end do

      ! loop through the input files and weights
      do ndx = 1, nfile
        filename = infile(ndx)
        write(*,*) 'Input: Opening ', filename
        open(2,FILE=filename(1:len_trim(filename)),STATUS='OLD')

        ! set histogram arrays so end is filled where not read
        do idx = 1, 12
          do jdx = 1, 16
            do kdx = 1, nbin
              freq(idx,jdx,kdx)=1.0
            end do
          end do
        end do

        !read everything in
        read (2,10) pound, wmoid, country, state, name
10      format (a2, a6,x, a2,x,a2,x,a20)

       !read station latitude, longitude, and elevation (m)
        read (2,*) lat_deg, lat_min, lat_NS, lon_deg, lon_min, lon_EW, elev
        latitude = lat_deg + lat_min/60.0
        if (lat_NS=='S') latitude = -1*latitude
        longitude = lon_deg + lon_min/60.0
        if (lon_EW=='W') longitude = -1*longitude

        ! store this point while converting lat lon degrees to radians
        trx(ndx) = sc * longitude
        try(ndx) = sc * latitude

        !read direction frequencies
        do idx=1,16
          read (2,*) (dirf(jdx,idx),jdx=1,12)
        end do
        !read calm frequencies
        read (2,*) (calm(idx),idx=1,12)

        !read all histogram info
        do idx = 1, 12
          do 22 jdx = 1, 16
            do kdx = 1, nbin
              read(2,49,EOR=22,ADVANCE='NO') freq(idx,jdx,kdx)
49            format(f6.3)
            end do
22        continue
        end do
        !max/min ratio
        read (2,*) (maxminr(idx),idx=1,12)
        !hour of the day when max. occurs
        read (2,*) (hrmax(idx),idx=1,12)

        ! close file for reading
        close(2)

        ! add read values to accumulators
        welev = welev + elev * weight(ndx)
        do idx = 1, 12
          wcalm(idx) = wcalm(idx) + calm(idx) * weight(ndx)
          wmaxminr(idx) = wmaxminr(idx) + maxminr(idx) * weight(ndx)
          whrmax(idx) = whrmax(idx) + hrmax(idx) * weight(ndx)
          ! round hour of maximum wind to nearest integer
          iwhrmax(idx) = int(whrmax(idx) + 0.4999999999)
          do jdx = 1, 16
            ! direction arrays
            wdirf(idx,jdx) = wdirf(idx,jdx) + dirf(idx,jdx) * weight(ndx)
            do kdx = 1, nbin
              ! histogram arrays
              wfreq(idx,jdx,kdx) = wfreq(idx,jdx,kdx) + freq(idx,jdx,kdx) * weight(ndx)
            end do
          end do
        end do

      end do

      ! Transform point lat/lon to Cartesian coordinates 
      call trans(nfile, try, trx, trx, try, trz)
      ! interpolate coordinates to point of interest
      px = 0.0
      py = 0.0
      pz = 0.0
      do ndx = 1, nfile
        px = px + trx(ndx) * weight(ndx)
        py = py + try(ndx) * weight(ndx)
        pz = pz + trz(ndx) * weight(ndx)
      end do
      ! Transform from cartesian back to lat/lon coordinates
      call scoord(px, py, pz, wlatitude, wlongitude, vnrm)
      ! convert back to degrees
      wlatitude = wlatitude/sc
      wlongitude = wlongitude/sc

      write(*,*) 'Output: Opening ', outfile
      open(1,FILE=outfile(1:len_trim(outfile)),STATUS='REPLACE')

      !write header line
      pound = '# '
      wmoid = '999999'
      country = 'US'
      state = 'XX'
      name = 'INTERPOLATED'
      write (1,110) pound, wmoid, country, state, name
110   format (a2, a6,x, a2,x,a2,x,a20)

      ! write latitude, longitude, and elevation (m)
      if (wlatitude .lt. 0.0) then 
        lat_NS = 'S'
        wlatitude = -1 * wlatitude
      else
        lat_NS = 'N'
      end if
      if (wlongitude .lt. 0.0) then 
        lon_EW = 'W'
        wlongitude = -1 * wlongitude
      else
        lon_EW = 'E'
      end if
      ! find degrees and minutes and round to nearest minute
      lat_deg = int(wlatitude + 0.5/60.0)
      lat_min = int((wlatitude - lat_deg) * 60.0)
      lon_deg = int(wlongitude + 0.5/60.0)
      lon_min = int((wlongitude - lon_deg) * 60.0)
      elev = int(welev + 0.5)
      !write (1,*) lat_deg,' ',lat_min,' ',lat_NS,' ',lon_deg,' ',lon_min,' ',lon_EW,' ',elev,' 20090101 20090101 XXX'
      write (1,*) lat_deg,lat_min,lat_NS,lon_deg,lon_min,lon_EW,elev,'20090101 20090101 XXX'

      ! write direction frequencies
      do idx=1,16
        write (1,120) (wdirf(jdx,idx),jdx=1,12)
120     format (12f5.1)
      end do
      !read calm frequencies
      write (1,120) (wcalm(idx),idx=1,12)

      !read all histogram info
      do idx = 1, 12
        do jdx = 1, 16
          do kdx = 1, nbin
            if (wfreq(idx,jdx,kdx) .lt. 1.0) then
              ibin = kdx
            else
              exit
            end if
          end do
          write(1,130) (wfreq(idx,jdx,kdx),kdx=1,ibin)
130       format(f5.3, 24f6.3)
        end do
      end do
      !max/min ratio
      write (1,120) (wmaxminr(idx),idx=1,12)
      !hour of the day when max. occurs
      write (1,140) (iwhrmax(idx),idx=1,12)
140   format (12i5)      


end