!
      module m_rdata_boxes
!
! Module to perform data thinning by selecting a single "best"
! observation for each geogaphical "box." It essentially saves 
! the complete header information required to simulate the 
! observation in an array named "box_info." An additional corresponding 
! array named "box_spot" is required in the case of AQUA AIRS and AMSU 
! observations to hold the "AQUASPOT" header information.  
!
! The thinning boxes defined here need not coincide with those employed
! in the data assimilation system.  In fact, it is likely better that 
! several of these here fall into one of the latter. In this way, 
! the data assimilation still performs its own data selection (but on 
! a reduced number of considered observations), yet the number of 
! simulated observation produced for the OSSE is much reduced compared to
! all available (with most then simply discarded by the observation 
! thinning algorithm employed by the data assimilation system).
!
! Initial code by Ronald Errico January 2008
!
!    History
!    2008-09-18 yang - modify to include msu
      use m_kinds
!
      implicit none      
!
      private
      public :: clean_rdata_boxes
      public :: fill_prof_info
      public :: find_box_index 
      public :: find_rad_subset
      public :: get_box_info 
      public :: get_box_spot_info 
      public :: put_box_info 
      public :: put_box_spot_info 
      public :: set_box_skip
      public :: set_rdata_boxes
!
      integer,                   save :: n_hdr
      integer,                   save :: n_hdr_sat
      integer,                   save :: n_hdr4
      integer,                   save :: n_types
      integer,                   save :: nlats
      integer,      allocatable, save :: nlonP(:)
      integer,      allocatable, save :: nlons(:)
      real(rkind1),              save :: dlat  
      real(rkind1), allocatable, save :: dlons(:) 
      real(rkind1), allocatable, save :: box_info(:,:,:)
      real(rkind3), allocatable, save :: box_spot(:,:)
!
      contains
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine clean_rdata_boxes 
      implicit none
!
      deallocate (nlonP,dlons,nlons,box_info)
      print *,'Deallocate: nlonP, dlons, box_info'
!
      if (n_hdr_sat /= 0) then
        deallocate (box_spot)
        print *,'Deallocate: box_spot'
      endif
!
      end subroutine clean_rdata_boxes
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine compute_hours (datetime,hours,lerror)
!
      implicit none 
      integer, parameter :: rkind8=8
      integer :: datetime(6)
      real(rkind8) :: hours
      logical :: lerror
!
      logical :: leapyear
      integer :: i
      integer :: iyear
      integer :: lmday
      integer :: mdays(12)
      integer :: nleaps
      real(rkind8), parameter :: zero=0._rkind8
      real(rkind8), parameter :: rsixty=1._rkind8/60._rkind8
      data mdays /31,28,31,30,31,30,31,31,30,31,30,31/
!
      lerror=.false.
      hours=zero  
!
      iyear=datetime(1)-2001
      if (iyear < 0  ) then 
        lerror=.true.
      else
        hours=hours+iyear*24.*365.
      endif 
!
      if (mod(iyear,4) == 3) then
        leapyear=.true.
      else
        leapyear=.false.
      endif
!
      if ( (datetime(2) < 1) .or. (datetime(2) > 12) ) then
        lerror=.true.
      else
        if (datetime(2) > 1) then
          do i=1,datetime(2)-1
            hours=hours+24.*mdays(i)          
          enddo
        endif
        nleaps=(iyear+1)/4
        if ((nleaps > 1) .and. leapyear .and. (datetime(2) ==1 )) then
          nleaps=nleaps-1
        endif
        hours=hours+nleaps*24.
      endif
!
      if ( (datetime(2) == 2) .and. leapyear) then
        lmday=29 
      else
        lmday=mdays(datetime(2))
      endif
      if ((datetime(3) < 1) .or. (datetime(3) > lmday) ) then
        lerror=.true.
      else
        hours=hours+24.*(datetime(3)-1.)
      endif
!
      if ( (datetime(4) < 0) .or. (datetime(4) > 24) ) then 
        lerror=.true.
      else
        hours=hours+datetime(4)
      endif
!
      if ( (datetime(5) < 0) .or. (datetime(5) > 60) ) then 
        lerror=.true.
      else
        hours=hours+rsixty*datetime(5)
      endif
!
      if ( (datetime(6) < 0) .or. (datetime(6) > 60) ) then 
        lerror=.true.
      else
        hours=hours+rsixty*rsixty*datetime(6)
      endif

!
      end subroutine compute_hours 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine fill_prof_info (idatetime,n_errors,ierrors,obs_info, &
                                 d_type,prof_info,lerror)
!              
      implicit none
!
      integer*8 :: idatetime
      integer :: n_errors
      integer :: ierrors(n_errors)
      real(rkind3) :: obs_info(n_hdr)
      real(rkind1) :: prof_info(n_hdr4)
      character(len=*) :: d_type
      logical :: lerror
!
      logical :: lerrors1, lerrors2
      integer :: datetime(4)
      integer :: dt6(6)
      integer :: i
      integer*8 :: id
      real(rkind3) :: dhours
      real(rkind3) :: hours_syn
      real(rkind3) :: hours_obs
      real(rkind3) :: clat
      real(rkind3) :: clon
      real(rkind3), parameter :: r90=90._rkind3
      real(rkind3), parameter :: r180=180._rkind3
      real(rkind3), parameter :: r360=360._rkind3
!
!
!  Copy to prof_info
      prof_info(1:n_hdr)=obs_info(1:n_hdr)
!
      if (obs_info(2) > 1.e9) then  ! missing data
        prof_info(n_hdr+4)=1.e10
        lerror=.true.
!
      else                          ! treatment of non-missing data
!
! unpack datetime
        datetime(1)=idatetime/1000000
        id=idatetime-datetime(1)*1000000
        datetime(2)=id/10000
        id=id-datetime(2)*10000
        datetime(3)=id/100
        datetime(4)=id-datetime(3)*100
!
! Compute hours from start of 2001 for requested date/time group      
        dt6(1:4)=datetime(1:4)
        dt6(5:6)=0
        call compute_hours (dt6,hours_syn,lerrors1)
!  
! Compute hours from start of 2001 for observation info
        do i=1,6 
          dt6(i)=nint(obs_info(i+3))
        enddo
        call compute_hours (dt6,hours_obs,lerrors2)
!
! Compute departure from synoptic time or flag as error
        if (lerrors1.or.lerrors2) then
          lerror=.true.
        else
          lerror=.false.  
          dhours=hours_obs-hours_syn
        endif
!
! Check if time range of observation is acceptable
        if (dhours > 3.) then
          ierrors(1)=ierrors(1)+1
          lerror=.true.
        endif
        if (dhours < -3.) then
          ierrors(2)=ierrors(2)+1
          lerror=.true.
        endif
        prof_info(n_hdr+4)=dhours     
!
! Check lat and lon 
        clon=obs_info(3)
        if ( (clon < -r180) .or. (clon > r360) ) then
          ierrors(3)=ierrors(3)+1
          lerror=.true.
        endif
!
        clat=obs_info(2)
        if ( (clat < -r90) .or. (clat > r90) ) then
          ierrors(4)=ierrors(4)+1
          lerror=.true.
        endif
! 
      endif  ! test if missinig data
!
      end subroutine fill_prof_info 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine find_box_index (lat,lon,index)
!
! find filtering-box index
! 
      implicit none
      real(rkind1) :: lat
      real(rkind1) :: lon
      integer      :: index
!
      integer :: nlat, nlon
      real(rkind1) :: lonx
!
      nlat=1+int((lat+90.)/dlat)
      nlat=min(nlat,nlats)
!
      if (lon < zero_k1 ) then
        lonx=360.+ lon
      else
        lonx=lon
      endif
      nlon=1+int(lonx/dlons(nlat))
      nlon=min(nlon,nlons(nlat))
!
      index=nlon+nlonP(nlat)
!
      end subroutine find_box_index 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine find_rad_subset (prof_ndim4,prof_info,d_type, &
                                  ierrors,lerror,id)
!
! Distinguish radiance observation subsets based on Satellite ID if not 
! AIRS, otherwise instrument ID if AIRS.  Return id=1 if the subset is
! not listed here as one of the acceptable ones.
!
! For AIRS, the file observational data is assumed to contain both the
! AIRS instrument and the AQUA-AMSUA data.
!
! For HIRS2, only one satellite is used; for HIRS3, only two.
! For MSU, only one NOAA-14 is used.
!
! Note that if this routine requires changes, likely also does the
! subroutine "read_error_table" used with the program "add_error"
!
      use m_kinds
      implicit none
      logical :: lerror
      integer :: id
      integer :: ierrors
      integer :: prof_ndim4
      real(rkind1) :: prof_info(prof_ndim4)
      character(len=*) :: d_type
!
      integer, parameter :: i_siid=1
      integer, parameter :: i_said=14
!
      id=0
      if (d_type == 'HIRS2') then
        if (prof_info(i_said) == 205) id=1  ! only NOAA 14
      elseif (d_type == 'HIRS3') then         
        if (prof_info(i_said) == 207) id=1  ! NOAA 16
        if (prof_info(i_said) == 208) id=2  ! NOAA 17
      elseif (d_type(1:4) == 'AMSU') then
        if (prof_info(i_said) == 205) id=1  ! NOAA 14
        if (prof_info(i_said) == 206) id=2  ! NOAA 15
        if (prof_info(i_said) == 207) id=3  ! NOAA 16
        if (prof_info(i_said) == 208) id=4  ! NOAA 17
        if (prof_info(i_said) == 209) id=5  ! NOAA 18
      elseif (d_type == '_MSU_') then
        if (prof_info(i_said) == 205) id=1  ! NOAA 14
      else ! AIRS
        if (prof_info(i_siid) == 420) id=1  ! AQUA AIRS IR
        if (prof_info(i_siid) == 570) id=2  ! AQUA AMSU-A
      endif
!
      if (id == 0) then
        lerror=.true.
        ierrors=ierrors+1
      endif
!
     if (id > n_types) then
       print *,' '
       print *,' Error found in subroutine find_rad_subset'
       print *,' id=',id,' greater than n_types=',n_types
       stop
     endif 
!
      end  subroutine find_rad_subset 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine get_box_info (index,itype,b_info)
!
! Retrieve information from thinning box
      implicit none
      integer :: index
      integer :: itype
      real(rkind1) :: b_info(n_hdr4)
      b_info(1:n_hdr4)=box_info(1:n_hdr4,index,itype)
      end subroutine get_box_info 
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine get_box_spot_info (index,b_info)
!
! Retrieve information from box holding additional "SPOT" information
      implicit none
      integer :: index
      real(rkind3) :: b_info(n_hdr_sat)
      b_info(1:n_hdr_sat)=box_spot(1:n_hdr_sat,index)
      end subroutine get_box_spot_info 
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine put_box_info (index,itype,lcldtop,pcldtop,b_info,l_box)
!
! Place (or over-write) information in a thinning box
!
      implicit none
      integer :: index     ! index of box number (location)
      integer :: itype     ! index of observation sub-type
      integer :: lcldtop   ! number of profile p-levles from top to cloud
      real(rkind1) :: pcldtop            ! pressure of cloud top
      real(rkind1) :: b_info(n_hdr4)     ! scalar info for obs spot
      logical :: l_box     ! .true. if box_info over written
      integer :: icol
!
! the following conditional tests whether
!  (1) this is the first obs falling in the thinning box, and therefore
!      this obs's scalar info should be saved, or
!  (2) this obs is in a clearer profile than any previous obs considered
!      for this thinning box (i.e., any clouds are located lower in the 
!      profile), and therefore this new obs should be saved in place of 
!      any earlier one.
!  (3) if the same number of levels are cloud-free, then retain the obs.
!      closer to the synoptic time.
!  (4) a value of lcldtop=1 serves as a flag to omit a record
!
      l_box=.false. 
      if (lcldtop > 1) then
        icol=n_hdr+1
        box_info(icol,index,itype)=box_info(icol,index,itype)+one_k3
        if ( (nint(box_info(icol,index,itype)) == 1) .or.               &
             (nint(box_info(icol+1,index,itype)) < lcldtop) .or.        &
             ( (nint(box_info(icol+1,index,itype)) == lcldtop) .and.    &
               (abs(box_info(icol+3,index,itype)) > abs(b_info(icol+3)) &
               ) ) ) then 
          box_info(icol+1,index,itype)=lcldtop
          box_info(icol+2,index,itype)=pcldtop
          box_info(1:n_hdr,index,itype)=b_info(1:n_hdr)
          box_info(icol+3:n_hdr4,index,itype)=b_info(icol+3:n_hdr4)
          l_box=.true.
        endif
      endif   ! test on lcldtop
!
!
      end subroutine put_box_info 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine put_box_spot_info (index,b_info)
!
! Place information in a box to save additional "SPOT" information
! For AQUA, this includes combined header info for AIRS, AMSUA, and HSB

      implicit none
      integer :: index
      real(rkind3) :: b_info(n_hdr_sat)
      box_spot(1:n_hdr_sat,index)=b_info(1:n_hdr_sat)
      end subroutine put_box_spot_info 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine set_box_skip (nboxes,ntypes,box_skip,nums)
!
!  Count the number of non-empty boxes (num)
!  For each individual subtype, shift the content of boxes so that the 
!   first num boxes are the non-empty ones and the remaining boxes empty
!  Determine indexes of those subtypes that have some non-empty boxes
!
      implicit none
      integer :: nboxes
      integer :: ntypes
      integer :: box_skip(ntypes)
      integer :: nums(ntypes) 
!
      integer :: i, j, j1
!
      if (ntypes /= n_types) then
        print *,' '
        print *,'Error in subroutine set_box_skip'
        print *,'ntypes=',ntypes,' differs from n_types=',n_types
        stop
      endif
!
! Count number of non-empty thinning boxes for each subtype
! Also, shift contents so the first nums(j) boxes are the non-empty ones
      nums(:)=0       
      do j=1,ntypes
        do i=1,nboxes
          if (nint(box_info(n_hdr4-3,i,j)) /= 0) then
            nums(j)=nums(j)+1
            if (i /= nums(j)) then
              box_info(:,nums(j),j)=box_info(:,i,j)
              box_info(:,i,j)=zero_k1
! Shift contents of sat header ("spot") info stored in thinning boxes
              if ((j == 1) .and. (n_hdr_sat > 0)) then 
                box_spot(1:n_hdr_sat,nums(j))=box_spot(1:n_hdr_sat,i)
                box_spot(1:n_hdr_sat,i)=zero_k1
              endif
            endif
          endif  ! test on whether box is not empty
        enddo    ! loop over boxes
      enddo      ! loop over subtypes
! 
! Determine indexes for subtypes with non-empty boxes  
      j=0
      box_skip(:)=0
      do i=1,n_types
        if (nums(i) /= 0) then 
          j=j+1
          box_skip(j)=i
        endif        
      enddo
!
      print *,' '
      print *,' Numbers of profiles to be considered for each subtype:'
      print ('(6i7)'),nums(1:ntypes)
      print *,' Indexes of detected subtypes:'
      print ('(6i7)'),box_skip(1:ntypes)
!
      end subroutine set_box_skip 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X 
!
      subroutine set_rdata_boxes (box_size,n_spot,ndim4,nboxes,ntypes, &
                                  n_spot2) 
!
      implicit none
!
      integer      :: nboxes
      integer      :: ntypes
      integer      :: n_spot   ! header size for instrument report
      integer      :: n_spot2  ! header size for "AQUASPOT" if needed
      integer      :: ndim4    ! size of array prof_info
      real(rkind1) :: box_size
!
      integer      :: n
      real(rkind1), parameter :: earth_r=6370.  ! radius of earth in km
      real(rkind1) :: pi, lat, circum
!
      pi=4.*atan(1.)
      nlats=int(pi*earth_r/box_size)
      dlat=180./nlats
      allocate (nlonP(nlats))
      allocate (nlons(nlats))
      allocate (dlons(nlats))
!      
      nboxes=0
      do n=1,nlats
        lat=-90.+(n-.5)*dlat
        nlonP(n)=nboxes
        circum=2.*pi*earth_r*cos(pi*lat/180.)
        nlons(n)=max(1,int(circum/box_size))
        dlons(n)=360./nlons(n)
        nboxes=nboxes+nlons(n)
      enddo
!
      n_hdr=n_spot      ! save this value
      n_hdr_sat=n_spot2 ! save this value
      n_hdr4=ndim4      ! save this value
      n_types=ntypes    ! save this value
!
      allocate (box_info(n_hdr4,nboxes,ntypes))
      box_info(1,:,:)=0.
!
      print *,' '
      print ('(a,i10,a)'),' Thinning boxes defined for ',nboxes,' boxes'
      print ('(a,f8.1,a,i4,f6.2,a,i3)'),' box_size=',box_size, &
                  ', nlats,dlat=',nlats,dlat,', ntypes=',ntypes 
!
! Create additional space for saving additional "SPOT" info 
! corresponding to box_info
      if (n_spot2 > 0) then
        allocate (box_spot(n_spot2,nboxes))
        print *,' Additional thinning box created for storing', &
                ' satellite spot info:'
        print ('(a,i3,a,i10)'),' n_spot2=',n_spot2,', nboxes=',nboxes 
      endif
      
!
      end subroutine set_rdata_boxes
!
!
      end module m_rdata_boxes