subroutine reorg_metar_cloud(cdata,nreal,ndata,cdata_all,maxobs,ngrid)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:  reorg_metar_cloud     define a closest METAR cloud observation for each grid point
!
!   PRGMMR: Ming Hu          ORG: GSD/AMB        DATE: 2009-09-21
!
! ABSTRACT: 
!
! PROGRAM HISTORY LOG:
!    2010-04-21  Hu  initial
!
!  input argument list:
!     cdata     - METRA cloud observation
!     nreal     - first dimension of cdata
!     ndata     - second dimension of cdata
!     maxobs    - maximum number of cdata_all
!     ndata     - number of type "obstype" observations retained for further processing
!
!  output argument list:
!     cdata_all - reorganized METAR observation
!     ngrid     - number of type metar cloud observations reorganized to the grid
!
! USAGE:
!   INPUT FILES: 
!   OUTPUT FILES:
!
! REMARKS:
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90 
!   MACHINE:  Linux cluster (WJET)
!
!$$$
!
!_____________________________________________________________________
!

  use kinds, only: r_kind,i_kind,r_double
  use gridmod, only: nlon,nlat
  use constants, only: one,half
  use rapidrefresh_cldsurf_mod, only: metar_impact_radius

  implicit none

! Declare passed variables
  integer(i_kind)                       ,intent(in) :: nreal   
  integer(i_kind)                       ,intent(in) :: ndata
  integer(i_kind)                       ,intent(in) :: maxobs
  real(r_kind),dimension(nreal,ndata)   ,intent(inout) :: cdata
  real(r_kind),dimension(nreal,maxobs)  ,intent(out):: cdata_all
  integer(i_kind)                       ,intent(out):: ngrid

! local variable
!
  integer(i_kind) :: ista_prev,ista_prev2,ista_save

  INTEGER(i_kind),allocatable :: first_sta(:,:)
  INTEGER(i_kind),allocatable :: next_sta(:)
  INTEGER(i_kind) ::    null_p
  PARAMETER ( null_p     = -1       )

  INTEGER(i_kind)  :: nsta_cld
  INTEGER(i_kind),allocatable :: sta_cld(:)


  integer(i_kind) :: isprd,isprd2,iout
  integer(i_kind) ::  aninc_cld_p
  integer(i_kind) :: i,j,k,i1,j1,ic,ista
  integer(i_kind) ::  ista_min
  real(r_kind) ::  min_dist, dist

  real(r_double) rstation_id
  character(8) :: cstation1,cc,ci
  equivalence(cstation1,rstation_id)
!
! 
  isprd=int(metar_impact_radius + half)
  if(isprd <= 0) return
  if(ndata <= 0) return
  ngrid = 0
  aninc_cld_p = isprd
  isprd2      = isprd * isprd
!
!  check duplicate observation and pick the one near analyis time
!
  DO i=1,ndata-1
    DO ic = i+1, ndata
      rstation_id=cdata(1,ic)
      cc=cstation1
      rstation_id=cdata(1,i)
      ci=cstation1
      if(trim(cc) == trim(ci) ) then
        if(abs(cdata(21,ic)) >= abs(cdata(21,i))) then   ! 21= observation time
          cdata(22,ic)=100                               ! 22= usage
        else
          cdata(22,i)=100
        endif
      endif
    ENDDO  !  ic
  ENDDO  !  i
!
!  find the nearest METAR station for each grid point
!
  allocate(first_sta(nlon,nlat))
  allocate(next_sta(ndata))
  allocate(sta_cld(ndata))
  first_sta= null_p
  next_sta = null_p
  sta_cld  = 0

! -- Then fill arrays based on where obs are
! -- NEXT_STA  - pointer to next ob found in grid volume
!         from the previous ob in that volume
!-----------------------------------------------------------
  do ista = 1,ndata
     if (cdata(22,ista) < 50) then
       if (cdata(6,ista) >= -one .and. cdata(6,ista) < 100.0_r_kind) then
! when cloud data are available 

          i = int(cdata(2,ista))
          j = int(cdata(3,ista))

          if (first_sta(i,j) == null_p) then
             first_sta(i,j) = ista
          else
             ista_prev = first_sta(i,j)
             do while (ista_prev /= null_p )
                ista_prev2= next_sta(ista_prev)
                ista_save = ista_prev
                ista_prev = ista_prev2
             enddo
             next_sta(ista_save) = ista
          end if
       endif
     endif
   enddo

   iout=0
   DO j = 1,nlat,aninc_cld_p
     DO i = 1,nlon,aninc_cld_p
       nsta_cld = 0

!sb -- Find all stations w/ cloud data within encompassing box
!mh   The box is decide by the influence radius of the analysis   
       do j1 = max(1,j-isprd),min(nlat,j+aninc_cld_p+isprd)
         do i1 = max(1,i-isprd),min(nlon,i+aninc_cld_p+isprd)
           if (first_sta(i1,j1) /= null_p) then
              ista_prev = first_sta(i1,j1)
              do while (ista_prev /= null_p )
                 nsta_cld = nsta_cld + 1
                 sta_cld(nsta_cld) = ista_prev            
                 ista_prev2= next_sta(ista_prev)
                 ista_prev = ista_prev2
              enddo
           end if
         enddo   ! j1
       enddo   ! i1

!
!sb - This is the big grid pt loop.  Walk thru each
!            individual grid points within 10x10 box.
       if(nsta_cld > 0 ) then
         do j1 = j,min(nlat,j+aninc_cld_p - 1)
           do i1 = i,min(nlon,i+aninc_cld_p - 1)

!sb - Find closest cloud station to grid point
             min_dist = 1.e10_r_kind
             do ic= 1,nsta_cld
                ista = sta_cld(ic)
                dist = (float(i1)-cdata(2,ista))*(float(i1)-cdata(2,ista))  &
                      +(float(j1)-cdata(3,ista))*(float(j1)-cdata(3,ista))
                if (dist < min_dist .and. dist < float(isprd2)) then
                    min_dist = dist
                    ista_min = ista
                end if
             end do  ! ic find the closest cloud station

             if (min_dist < 1.e9_r_kind) then
                iout = iout + 1
                if(iout > maxobs) then
                  write(6,*)'reorg_metar_cloud:  ***Error*** ndata > maxobs '
                  call stop2(50)
                end if
                do k=1,nreal
                   cdata_all(k,iout) = cdata(k,ista_min)
                enddo
                cdata_all(2,iout) = float(i1)        ! grid index i
                cdata_all(3,iout) = float(j1)        ! grid index j
                cdata_all(nreal,iout)=sqrt(min_dist) ! distance from station
             endif
           enddo   ! j1
         enddo   ! i1
       endif ! nsta_cld > 0
     enddo   ! j
   enddo   ! i

   ngrid = iout

   deallocate(first_sta)
   deallocate(sta_cld)
   deallocate(next_sta)

end subroutine reorg_metar_cloud