subroutine da_create_bins(ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
   lat_min, lat_max, binwidth_lat, hgt_min, hgt_max, binwidth_hgt, latitude, height)

   !----------------------------------------------------------------------------
   !
   ! Purpose: To create the bins for calculation of statistics.
   !
   ! Input:
   ! ni, nj, nk   Dimensions
   ! bin_type     0: No binning; 
   !              1: bins for X-direction mean; 
   !              2: bins for each of binwidth_lat/binwidth_hgt.  
   !              3: bins for each of binwidth_lat/nk.  
   !              4: num_hor_bins horizontal bins /nk.  
   !              5: Average over all horizontal points (nk bins for 3D fields)
   !              6: Average over all points (only 1 bin).
   ! Optional for bin_type = 2:
   ! lat_min, lat_max Minimum/maximum latitude ranges for bin_type = 2
   ! binwidth_lat interval between bins for latitude in degree for bin_type = 2
   ! binwidth_hgt interval between bins for height in meter for bin_type = 2
   ! num_bins_hgt the number of height bins for bin_type = 2
   ! latitude     3d field latitude in degree for bin_type = 2
   ! height       3d field height in meter for bin_type = 2
   !
   ! Output:
   ! num_bins,num_bins2d ---- the number of bins for 3d and 2d fields
   ! bin     ,     bin2d ---- Assigned bin to a gridpoints for 3d and 2d fields
   !
   !----------------------------------------------------------------------------

   implicit none

   integer,        intent(in)  :: ni, nj, nk          ! Dimensions read in.
   integer,        intent(in)  :: bin_type            ! Type of bin to average over
   integer,        intent(out) :: num_bins            ! Number of bins.
   integer,        intent(out) :: num_bins2d          ! Number of bins for 2D fields
   integer,        intent(out) :: bin(1:ni,1:nj,1:nk) ! Bin at each 3D point
   integer,        intent(out) :: bin2d(1:ni,1:nj)    ! Bin at each 2D point

   real, optional, intent(in)  :: lat_min, lat_max    ! Used if bin_type = 2 (deg)
   real, optional, intent(in)  :: binwidth_lat        ! Used if bin_type = 2 (deg)
   real, optional, intent(in)  :: hgt_min, hgt_max    ! Used if bin_type = 2 (deg)
   real, optional, intent(in)  :: binwidth_hgt        ! Used if bin_type = 2 (m).
   real, optional, intent(in)  :: latitude(1:ni,1:nj) ! Latitude (degrees).
   real, optional, intent(in)  :: height(1:ni,1:nj,1:nk)     ! Height (m).

   integer           :: b, i, j, k                 ! Loop counters.
   integer           :: count                      ! Counter
   integer           :: num_bins_lat               ! Used if bin_type = 2.
   integer           :: num_bins_hgt               ! Used if bin_type = 2.
   integer           :: bin_lat                    ! Latitude bin.
   integer           :: bin_hgt                    ! Height bin.
   integer           :: num_bins_i, num_bins_j     ! Used if bin_type = 4.
   integer           :: nii, njj                   ! Used if bin_type = 4.
   integer           :: bin_i(1:ni), bin_j(1:nj)   ! Used if bin_type = 4.
   real, allocatable :: binstart_lat(:)            ! Used if bin_type = 2 (deg)
   real, allocatable :: binstart_hgt(:)            ! Used if bin_type = 2 (deg)

   if (trace_use) call da_trace_entry("da_create_bins")

   if (bin_type == 0) then         ! No averaging in space

      num_bins = nk * nj * ni
      num_bins2d = nj * ni    ! Equals number of horizontal points.

      count = 1
      do k = 1, nk
         do j = 1, nj
            do i = 1, ni
               bin(i,j,k) = count
               count = count + 1
            end do
         end do
      end do
      bin2d(:,:) = bin(:,:,1)

   else if (bin_type == 1) then    ! Average over x-direction.

      num_bins = nj * nk
      num_bins2d = nj

      count = 1
      do k = 1, nk
         do j = 1, nj
            bin(1:ni,j,k) = count
            count = count + 1
         end do
      end do
      bin2d(:,:) = bin(:,:,1)

   else if (bin_type == 2) then    ! Global latitude/height bins:

      ! Setup latitude bins:
      write(unit=stdout,fmt='(/a,f12.5)')'   Minimum latitude = ', lat_min
      write(unit=stdout,fmt='(a,f12.5)')'    Maximum latitude = ', lat_max
      write(unit=stdout,fmt='(a,f12.5)') &
         '    Latitude bin width = ', binwidth_lat
      num_bins_lat = nint((lat_max - lat_min) / binwidth_lat)
      write(unit=stdout,fmt='(a,i8)') &
         '    Number of latitude bins = ', num_bins_lat
   
      allocate(binstart_lat(1:num_bins_lat))
      do b = 1, num_bins_lat ! Assume south to north (as in WRF).
         binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
      end do

      ! Setup height bins:
      write(unit=stdout,fmt='(/a,f12.5)')'    Minimum height = ', hgt_min
      write(unit=stdout,fmt='(a,f12.5)')'    Maximum height = ', hgt_max
      write(unit=stdout,fmt='(a,f12.5)')'    Height bin width = ', binwidth_hgt
      num_bins_hgt = nint((hgt_max - hgt_min) / binwidth_hgt)
      write(unit=stdout,fmt='(a,i8)') &
         '    Number of height bins = ', num_bins_hgt

      allocate(binstart_hgt(1:num_bins_hgt))
      do b = 1, num_bins_hgt
         binstart_hgt(b) = hgt_min + real(b-1) * binwidth_hgt
      end do

      num_bins = num_bins_lat * num_bins_hgt
      num_bins2d = num_bins_lat

      ! Select height bins:
      do j = 1, nj
         do i = 1, ni
            do k = 1, nk
               if (height(i,j,k) < binstart_hgt(1)) then 
                  bin_hgt = 1 ! In first bin.
               else if (height(i,j,k) >= binstart_hgt(num_bins_hgt)) then
                  bin_hgt = num_bins_hgt ! In final bin.
               else 
                  do b = 1, num_bins_hgt-1
                     if (height(i,j,k) >= binstart_hgt(b) .and. &
                          height(i,j,k) <  binstart_hgt(b+1)) then
                        bin_hgt = b
                        exit
                     end if
                  end do
               end if

               ! Select latitude bin that point falls in:
               if (k == 1) then
                  do b = 1, num_bins_lat-1
                     if (latitude(i,j) >= binstart_lat(b) .and. &
                        latitude(i,j) < binstart_lat(b+1)) then
                        bin_lat = b
                        exit
                     end if
                  end do
                  if (latitude(i,j) >= binstart_lat(num_bins_lat)) then
                     ! In final bin.
                     bin_lat = num_bins_lat
                  end if
                  bin2d(i,j) = bin_lat
               end if
               bin(i,j,k) = bin_lat + num_bins_lat * (bin_hgt - 1)
            end do
         end do
      end do

      deallocate(binstart_lat)
      deallocate(binstart_hgt)

   else if (bin_type == 3) then    ! Latitude/nk bins:

      ! Setup latitude bins:
      write(unit=stdout,fmt='(/a,f12.5)')'   Minimum latitude = ', lat_min
      write(unit=stdout,fmt='(a,f12.5)')'    Maximum latitude = ', lat_max
      write(unit=stdout,fmt='(a,f12.5)')'    Latitude bin width = ',binwidth_lat
      num_bins_lat = nint((lat_max - lat_min) / binwidth_lat)
      write(unit=stdout,fmt='(a,i8)') &
         '    Number of latitude bins = ', num_bins_lat
   
      allocate(binstart_lat(1:num_bins_lat))
      do b = 1, num_bins_lat ! Assume south to north (as in WRF).
         binstart_lat(b) = lat_min + real(b-1) * binwidth_lat
      end do

      num_bins = num_bins_lat * nk
      num_bins2d = num_bins_lat

      ! Select bins:
      do j = 1, nj
         do i = 1, ni
            do k = 1, nk
               ! Select latitude bin that point falls in:
               if (k == 1) then
                  do b = 1, num_bins_lat-1
                     if (latitude(i,j) >= binstart_lat(b) .and. &
                        latitude(i,j) < binstart_lat(b+1)) then
                        bin_lat = b
                        exit
                     end if
                  end do
                  if (latitude(i,j) >= binstart_lat(num_bins_lat)) then
                     ! In final bin.
                     bin_lat = num_bins_lat
                  end if
                  bin2d(i,j) = bin_lat
               end if
               bin(i,j,k) = bin_lat + num_bins_lat * (k - 1)
            end do
         end do
      end do

      deallocate(binstart_lat)

   else if (bin_type == 4) then    ! binwidth_lat/nk bins:

      ! Setup horizontal bins:
      write(unit=stdout,fmt='(/a,f12.5)') &
         '   Number of grid-cells to average over = ', binwidth_lat
      ! use binwidth_lat, but actually an integer number of points.
 
      num_bins_j = int(real(nj) / real(binwidth_lat))
      njj = int(binwidth_lat) * num_bins_j
      do j = 1, njj
         bin_j(j) = 1 + int(real(j-1) / binwidth_lat)
      end do
      if (nj > njj) bin_j(njj+1:nj) = bin_j(njj)

      num_bins_i = int(real(ni) / real(binwidth_lat))
      nii = int(binwidth_lat) * num_bins_i
      do i = 1, nii
         bin_i(i) = 1 + int(real(i-1) / binwidth_lat)
      end do
      if (ni > nii) bin_i(nii+1:ni) = bin_i(nii)

      num_bins2d = num_bins_i * num_bins_j
      num_bins = num_bins2d * nk

      do j = 1, nj
         do i = 1, ni
            bin2d(i,j) = bin_i(i) + (bin_j(j) - 1) * num_bins_i
            do k = 1, nk
               bin(i,j,k) = bin2d(i,j) + (k - 1) * num_bins2d
            end do
         end do
      end do

   else if (bin_type == 5) then    ! Average over all horizontal points.

      num_bins = nk
      num_bins2d = 1

      do k = 1, nk
         bin(:,:,k) = k
      end do
      bin2d(:,:) = 1

   else if (bin_type == 6) then    ! Average over all points.

      num_bins = 1
      num_bins2d = 1
      bin(:,:,:) = 1
      bin2d(:,:) = 1
   end if

   if (trace_use) call da_trace_exit("da_create_bins")

end subroutine da_create_bins