psmile_global_search_nnx_dble.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_global_search_nnx_dble
!
! !INTERFACE:

      subroutine psmile_global_search_nnx_dble (comp_info, search,        &
                        var_id, tgt_coords_x, tgt_coords_y, tgt_coords_z, &
                        neighbors_3d, nloc, num_neigh, nb_extra,          &
                        extra_search, send_index,                         &
                        ierror)
!
! !USES:
!
      use PRISM
!
      use psmile_grid, only : psmile_transform_extent_cyclic, &
                              psmile_transform_extent_back, &
                              max_num_trans_parts, &
                              code_no_trans, &
                              common_grid_range
!
      use PSMILe, dummy_interface => PSMILe_global_search_nnx_dble

      Implicit none
!
! !INPUT PARAMETERS:

      Type (Enddef_comp), Intent (In) :: comp_info
!
!     Info on the component in which the donor cells
!     should be searched.
!
      Type (Enddef_search), Intent (InOut) :: search

!     Info's on coordinates to be searched

      Integer,            Intent (In) :: var_id

!     Handle to the grid function

      Integer,            Intent (In) :: send_index

!     Send index of coupler send for the method and coordinates

      Integer,            Intent (In) :: nloc
!
!     Total  Number of locations to be transferred
!
      Double Precision,   Intent (In) :: tgt_coords_x (nloc)
      Double Precision,   Intent (In) :: tgt_coords_y (nloc)
      Double Precision,   Intent (In) :: tgt_coords_z (nloc)

!     Coordinates of the target points (which were searched and found)
!
      Integer,            Intent (In) :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.
!     besser n_corners

      Integer,            Intent (In) :: nb_extra
!
!     Number of nearest neigbhbours to be searched per extra point
!
! !INPUT/OUTPUT PARAMETERS:
!
      Integer,            Intent (InOut) :: neighbors_3d (ndim_3d, nloc, 
                                                          num_neigh)

!     Indices of neighbor locations (interpolation bases)
!
      Type (Extra_search_info), Intent (InOut) :: extra_search
!
!     Number of locations where
!       (*) required mask values were not true

! !OUTPUT PARAMETERS:
!
      Integer,           Intent (Out) :: ierror

!     Returns the error code of PSMILe_global_search_nnx_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS
!
!  DBLE_EARTH = Earth radius in meter
!
!  LEN_ITEM = Number of pieces of information
!             in integer buffer "ibuf"
!             1: index in indices/dist_dble/send_mask
!  INLP_LEN = Inner loop length for optimization
!  lon   = Index of Longitudes in arrays "boxes"
!  lat   = Index of Latitudes  in arrays "boxes"
!  vrt   = Index of height     in arrays "boxes"
!
!
      Double Precision, Parameter     :: dble_earth = 6400000.0d0
!
      Integer, Parameter              :: len_item = 1
!
      Integer, Parameter              :: inlp_len = 256
!
      Integer, Parameter              :: lon = 1
      Integer, Parameter              :: lat = 2
      Integer, Parameter              :: vrt = 3
!
!
! !LOCAL VARIABLES
!
!     ... Method pointer
!
      Integer                         :: comp_id, grid_id
      Integer                         :: global_grid_id
      Integer, Pointer                :: indices (:)
!
!     ... Info's on actual component
!
      Integer                         :: igrid, rank
      Integer                         :: n_total
      Integer                         :: grid_type
      Real (PSMILe_float_kind), Pointer :: extents (:, :, :)
!
!     ... lookup table
!
      Integer, allocatable            :: igrid_to_dest_comp(:)
!
!     ... loop variables
!
      Integer                         :: j, n
      Integer                         :: ibeg, iend
!
!     ... for data to other processes
!
      Logical                         :: mask_changed
      Integer                         :: n_extra
!
      Logical, Allocatable            :: send_mask (:)
      Double Precision, Allocatable   :: boxes (:, :, :)
      Real (PSMILe_float_kind)        :: range_box (2, ndim_3d)
      Real (PSMILe_float_kind)        :: dist_max (2)
      Double Precision                :: dble_huge
!
!     ... For locations controlled
!
!     Double Precision, Pointer       :: sin_search (:, :)
!     Double Precision, Pointer       :: cos_search (:, :)
!     Double Precision, Pointer       :: z_search (:)
      Double Precision, Pointer       :: dist_dble (:, :)
      Double Precision                :: r_earth2, r_lat
      Double Precision                :: delta, sin_d, cos_lat
!
!     ... transformed faces
!
      Integer                         :: n_trans, n_int
!
      Integer                         :: tr_codes (max_num_trans_parts)
      Integer                         :: found (max_num_trans_parts)
      Real (PSMILe_float_kind)        :: dinter_trans (2, ndim_3d)
      Real (PSMILe_float_kind)        :: transformed (2, ndim_3d, 
                                                      max_num_trans_parts), 
                                         dinter (2, ndim_3d, max_num_trans_parts)
!
      Double Precision, Save          :: period2 (ndim_3d)
!
!     ... for points found
!
      Integer                         :: n_found, n_liste
      Integer                         :: ip_dist
!
!     ... for send to partners
!
      Type (Send_information), Pointer :: send_info
      Integer                         :: extra_msg (nd_msgextra)
!
      Integer                         :: i, ip, ipi, ireq
      Integer                         :: len_ibuf, ndibuf, n_send
      Integer                         :: len_rbuf, len_rtem, ndrbuf
      Integer, Allocatable            :: ibuf (:)
      Double Precision, Pointer       :: buf (:)
!
!     ... for recv from partners
!
      Integer                         :: answer (nd_msgextra)
      Integer, Allocatable            :: selected (:, :)
      Type (Select_search_info), Pointer  :: sel_info (:)
!
!     ... for communication
!
      Integer                         :: dest, index, dest_comp, sender
      Integer                         :: nrecv, n_recv_req
      Integer                         :: recv_req
      Integer                         :: save_lreq (2:3)
      Integer                         :: status (MPI_STATUS_SIZE)
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 3
      Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_global_search_nnx_dble" performs a global 
! nearest neighbour search for the interpolation points 
! which require an extra seach.
!
! TODO: Save buffer indices in indices/send_mask/dist_dble (see len_item)
!       locally instead of passing them !!!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 09.10.06    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_global_search_nnx_dble.F90 2935 2011-02-03 07:52:41Z redler $
!  $Author: redler $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_global_search_nnx_dble.F90 2935 2011-02-03 07:52:41Z redler $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), var_id, search%msg_intersections%first_tgt_method_id, search%sender
#endif /* VERBOSE */

      call psmile_flushstd
!
      ierror = 0
!
      ndrbuf = 0
      ndibuf = 0
      n_recv_req = 0
!
      r_earth2 = 1.0d0 / (dble_earth * 2.0)
      r_lat = 1.0d0 / (dble_earth * dble_deg2rad)
!
      period2 = (common_grid_range(2,:) - common_grid_range(1,:)) / 2.0d0
!
      grid_id        = Methods(Fields(var_id)%method_id)%grid_id
      global_grid_id = Grids(grid_id)%global_grid_id
      grid_type      = Grids(grid_id)%grid_type
!
      comp_id = comp_info%comp_id
!
! Uebergeben oder ein-dimensional
!
      send_info => Methods(Fields(var_id)%method_id)%send_infos_coupler(send_index)
!
      indices    => extra_search%indices
      dist_dble  => extra_search%dist_dble
!
      n_extra = extra_search%n_extra
!
#ifdef PRISM_ASSERTION
#ifdef HUHU
      sin_search => extra_search%sin_search_dble
      cos_search => extra_search%cos_search_dble 
        z_search => extra_search%z_search_dble
      if ( .not. Associated (extra_search%dist_dble)       .or. &
           .not. Associated (extra_search%cos_search_dble) .or. &
           .not. Associated (extra_search%sin_search_dble) .or. &
           .not. Associated (extra_search%z_search_dble)   ) then
#else
      if ( .not. Associated (extra_search%dist_dble) ) then
#endif
!
         call psmile_assert (__FILE__, __LINE__, &
                             "arrays should be allocated and set")
      endif
#endif
!
!  Rank of actual process in Component communicator
!
      rank    = Comps(comp_id)%rank
!
! n_total   = Total number of grids/extents       in array "all_extents"
! igrid_end = Last  grid/extent of actual process in array "all_extents"
! igrid     = Index of actual grid in array "all_extents"
!
      extents => comp_info%all_extents
!
      n_total   = SUM (comp_info%Number_of_grids_vector(:))
!
      igrid = search%msg_intersections%first_src_all_extents_grid_id
!
! Can man die offset und partition information ausnutzen ?
!
!===> Allocate work vectors
!
! boxes(i)     = Box around ith extra point with interior radius
!                "dist_dble (i)"
! send_mask(i) = Does the i-th extra point to be transferred 
!                to the corresponding process ?
!
      Allocate (boxes (2, ndim_3d, n_extra), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = n_extra * (2 * ndim_3d)

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'boxes', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      Allocate (send_mask (n_extra), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = n_extra

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'send_mask', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      mask_changed = .true.
!
!===> Compute range (bounding box) for global nearest neighbour
!     search of extra points
!
!     ?dist_dble = huge ? neighbors_3d (:, i, nb_extra) auswerten ?
!     !!! TODO: Diese Punkte loesen eine Suche in allen Bloecken aus !
!               Besser erst lokal in der Nachbarschaft suchen !
!
!==> Get maximal extension for 3rd direction (z-direction)
!    ? einmal berechnen ?
!
     dble_huge = huge (dist_dble(1,1))
     dist_max (1) = extents (1, vrt, igrid)
     dist_max (2) = extents (2, vrt, igrid)
!
#ifdef DEBUGX
     print *, "dist_max", dist_max
#endif
!
      do n = 1, n_total
        if (global_grid_id == comp_info%all_extent_infos(4,n)) then
           dist_max (1) = min (dist_max(1), extents(1, vrt, n))
           dist_max (2) = max (dist_max(2), extents(2, vrt, n))
        endif
      end do
#ifdef DEBUGX
     print *, "dist_max", dist_max
#endif
!
      do ibeg = 1, n_extra, inlp_len
         iend = min (n_extra, ibeg + inlp_len - 1)
!
!        Construct boxes of interest in lon direction
!        Delta lon = 2 * asin (sin (dist_dble/(2*dble_earth)) / cos (lat))
!                    / dble_deg2rad
!cdir vector
         do i = ibeg, iend
            if (dist_dble(i, nb_extra) * r_earth2 >= dble_pih) then
               sin_d = 1.0d0
            else
               sin_d = abs(sin(dist_dble(i, nb_extra) * r_earth2))
            endif

            cos_lat = cos (tgt_coords_y(indices(i)) * dble_deg2rad)

            if (sin_d >= cos_lat) then
               delta = period2(lon)
            else
               delta = 2.0d0 * asin (sin_d / cos_lat) / dble_deg2rad
               delta = min (period2(lon), delta)
            endif

            boxes (1, lon, i) = tgt_coords_x(indices(i)) - delta
            boxes (2, lon, i) = tgt_coords_x(indices(i)) + delta
         end do

!cdir vector
         do i = ibeg, iend
            if (boxes (2,lon,i) - boxes (1,lon,i) > period2(lon)*2) then
               boxes (1,lon,i) = -period2(lon)
               boxes (2,lon,i) =  period2(lon)
            else if (boxes (1,lon,i) < -period2(lon)*2) then
               boxes (1,lon,i) = boxes (1,lon,i) + period2(lon)*2
               boxes (2,lon,i) = boxes (2,lon,i) + period2(lon)*2
            else if (boxes (2,lon,i) >  period2(lon)*2) then
               boxes (1,lon,i) = boxes (1,lon,i) - period2(lon)*2
               boxes (2,lon,i) = boxes (2,lon,i) - period2(lon)*2
            endif
         end do ! i

!        Construct boxes of interest in lat direction
!        Delta lat = dist_dble / (dble_earth * dble_deg2rad)

!cdir vector
         do i = ibeg, iend
            boxes (1, lat, i) = max (tgt_coords_y(indices(i)) -                        &
                                     min (period2(lat), dist_dble(i, nb_extra)*r_lat), &
                                     common_grid_range(1,lat))
            boxes (2, lat, i) = min (tgt_coords_y(indices(i)) +                        &
                                     min (period2(lat), dist_dble(i, nb_extra)*r_lat), &
                                     common_grid_range(2,lat))
         end do
      end do ! j
!
      range_box (1, lon) = minval (boxes (1, lon, :))
      range_box (2, lon) = maxval (boxes (2, lon, :))
!
      range_box (1, lat) = minval (boxes (1, lat, :))
      range_box (2, lat) = maxval (boxes (2, lat, :))
!
!        Construct boxes of interest in vertical direction
!cdir vector
      do i = 1, n_extra
         if (dist_dble(i, nb_extra) == dble_huge) then
            boxes (1, vrt, i) = dist_max (1)
            boxes (2, vrt, i) = dist_max (2)
         else
            boxes (1, vrt, i) = tgt_coords_z(indices(i)) -             &
                              dist_dble(i, nb_extra)
            boxes (2, vrt, i) = tgt_coords_z(indices(i)) +             &
                              dist_dble(i, nb_extra)
#ifdef DEBUGX
     print *, "dist_dble", i, indices(i), dist_dble(i, nb_extra)
#endif
         endif
      end do
!
      range_box (1, vrt) = minval (boxes (1, vrt, :))
      range_box (2, vrt) = maxval (boxes (2, vrt, :))
!
#ifdef PRISM_ASSERTION
      if (range_box (2,lon) - range_box (1,lon) > period2(lon)*4) then
          print *, 'range in lon direction', range_box (1:2,lon)
         call psmile_assert ( __FILE__, __LINE__, &
                             'range_box too large in lon direction')
      endif
!
      if (range_box (2,lat) - range_box (1,lat) > period2(lat)*4) then
          print *, 'range in lat direction', range_box (1:2,lat)
         call psmile_assert ( __FILE__, __LINE__, &
                             'range_box too large in lat direction')
      endif
#endif
!
!===> Transform range "range_box" into cyclic ranges
!
      call psmile_transform_extent_cyclic (grid_type,           &
                  range_box, transformed, tr_codes, n_trans, ierror)
      if (ierror > 0) return
!
!     create lookup table for translating igrid into a icomp
!
      Allocate (igrid_to_dest_comp(SUM(comp_info%Number_of_grids_vector)), stat = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = SUM(comp_info%Number_of_grids_vector)

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'igrid_to_dest_comp', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      igrid = 0
      do dest_comp = 1, comp_info%size
         do j = 1, comp_info%Number_of_grids_vector(dest_comp)
            igrid = igrid + 1
            igrid_to_dest_comp(igrid) = dest_comp
         enddo
      enddo
!
!===> Control extents of other domains in order to find partner
!
! igrid = current grid
!
!     loop over all grids/blocks (local and remote ones)
      do igrid = 1, SUM(comp_info%Number_of_grids_vector)

         dest_comp = igrid_to_dest_comp(igrid)

         if (mask_changed) then
            mask_changed  = .false.
            send_mask (:) = .false.
         endif
!
! TODO:
! Das ist mir zu teuer; erst mit den extents in einer einfachen
! schleife potentiell partner finden; siehe found (:)
! in psmile_find_intersect
!
!        if current block is not part of the same global grid
         if (global_grid_id /= comp_info%all_extent_infos(4,igrid)) cycle
!
!===> ... Is this my info
!         (note ranks start with 0, ...)
!
         if (rank == dest_comp-1 .and. &
               comp_info%all_extent_infos(1,igrid) == grid_id) cycle
!
!===> ... Does the "n"-th grid and "range_box" coincide ?
!        Rausziehen und vektorisieren ?
!
!  n_int = Number of significant intersections
!
!cdir vector
         do i = 1, n_trans
            dinter (1, :, i) = max (transformed(1,:,i), &
                                    extents(1,:,igrid))
            dinter (2, :, i) = min (transformed(2,:,i), &
                                    extents(2,:,igrid))
         end do
!
         n_int = 0
!cdir vector
         do i = 1, n_trans
            if (minval(dinter (2,:,i) - dinter (1,:,i)) >= 0.0d0) then
               n_int = n_int + 1
               found (n_int) = i
            endif
         end do
!
         if (n_int == 0) cycle
!
!===> ... Transform intersections back into original range
!
         do i = 1, n_int
            if (tr_codes(found(i)) /= code_no_trans) then
               call psmile_transform_extent_back (tr_codes(found(i)), &
                  dinter(:, :, found(i)), dinter_trans, 1, ierror)
               if ( ierror /= 0 ) return

               dinter (:, :, found(i)) = dinter_trans
            endif
         end do
!
!===> ... Get the extra points whose radius has an intersection
!         with the extent of the "n"-th grid
!
! send_mask(i) = Does the i-th extra point to be transferred 
!                to the corresponding process ?
!???  What about cyclic boundaries ????
!???  use global offsets here !
!
! Bedingung: 1. Punkt is in dinter drin
!            2. vereinfachte distanz kleiner als dist_dble
!
         do i = 1, n_int
#ifdef DEBUGX
print *, 'extent controlled', i, n_int, dinter(:,:,found(i))
#endif
!cdir vector
            do j = 1, n_extra
               send_mask (j) = send_mask (j) .or. &
                     (boxes(2,lon,j) >= dinter(1,lon,found(i)) .and. &
                      boxes(1,lon,j) <= dinter(2,lon,found(i)) .and. &
                      boxes(2,lat,j) >= dinter(1,lat,found(i)) .and. &
                      boxes(1,lat,j) <= dinter(2,lat,found(i)) .and. &
                      boxes(2,vrt,j) >= dinter(1,vrt,found(i)) .and. &
                      boxes(1,vrt,j) <= dinter(2,vrt,found(i)))
#ifdef DEBUGX
print *, j, boxes(:,:,j)
#endif
            end do ! j
         end do ! i
!
!===> ... Get number of points to be sent
!
         n_send = COUNT(send_mask)
         if (n_send == 0) cycle
!
         mask_changed = .true.
!
!===> ... A grid of process "dest_comp" is a potential partner
!
         dest = comp_info%psmile_ranks(dest_comp)
!
         n_recv_req = n_recv_req + 1
!
#ifdef VERBOSE
         print 9970, trim(ch_id), rank, dest_comp, dest, n_send

         call psmile_flushstd
#endif /* VERBOSE */
!
!===> ... Allocate send buffers
!
!  cf. subroutine PSMILe_Store_faces_3d_reg_dble
!
!  len_rtem = Number of pieces of information
!             per point in Double Precision buffer "rbuf"
!
!  Double Precision vector buf (1:ndrbuf):
!  buf (:) = Coordinates of the point to be searched
!            ndim_3d data items
!  buf (:) = Maximal distance
!            1 data item
!
!  Integer buffer ibuf (1:ndibuf):
!
!  ibuf (:) = Index in indices/dist_dble/send_mask
!
!        len_rtem = ndim_3d + nb_extra ! coords + alle distancen
         len_rtem = ndim_3d + 1        ! coords + groesste Distance
!
         len_rbuf = n_send * len_rtem
         len_ibuf = n_send * len_item
!
         if (len_rbuf > ndrbuf) then
            if (ndrbuf > 0) then
               Deallocate (buf)
            endif
!
            ndrbuf = len_rbuf
            Allocate (buf(ndrbuf), STAT = ierror)

            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = ndrbuf

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'buf', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
         endif
!
         if (len_ibuf > ndibuf) then
            if (ndibuf > 0) then
               Deallocate (ibuf)
            endif
!
            ndibuf = len_ibuf
            Allocate (ibuf(ndibuf), STAT = ierror)

            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = ndibuf

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'ibuf', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
         endif
!
!===> Store data in send buffers
!     ??? compact machen !?!?
!         Mehrfache Zellen Info's in ibuf und corner volumen
!         eliminieren ?
!
#ifdef USE_PACK
         ibuf (1:n_send*len_item:len_item) = &
            PACK ((/ (i, i=1,n_extra) /), send_mask)
#else
         ipi = 0

!cdir vector
         do i = 1, n_extra
            if (send_mask (i)) then
               ibuf (ipi+1) = i
               ipi = ipi + len_item
            endif
         end do
#endif /* USE_PACK */
!
!===> ... store coordinates and distance a single vectors 
!         (to be consistent with storage of PRISM external coordinates)
!         ? ist das wiklich besser, als die per Punkt zu speichern.
!           Grund fuer die Feldweise Speicherung ist, die
!           moegliche Wiederverwendung von Routinen.
!
! ? besser direkt cos_search, sin_search, z versenden ?
! Das waeren 5 satt 3 items;
! Neu rechnen schneller als senden ?
!              buf (ip+1) = cos_search (i, 1)
!              buf (ip+2) = cos_search (i, 2)
!              buf (ip+3) = sin_search (i, 1)
!              buf (ip+4) = sin_search (i, 2)
!              buf (ip+5) = z_search (i) = tgt_coords_z (indices(i))
!
#ifdef ALL_ITEM_PER_POINT
         ip = 0
!cdir vector
         do i = 1, n_extra
            if (send_mask (i)) then
!
               buf (ip+1) = tgt_coords_x (indices(i))
               buf (ip+2) = tgt_coords_y (indices(i))
               buf (ip+3) = tgt_coords_z (indices(i))
!
               buf (ip+4) = dist_dble (i, nb_extra)
               ip = ip + len_rtem
            endif
         end do
#else
         ip = 0
!cdir vector
         do i = 1, n_extra
            if (send_mask (i)) then
               ip = ip + 1
               buf (ip           ) = tgt_coords_x (indices(i))
               buf (ip+ n_send   ) = tgt_coords_y (indices(i))
               buf (ip+(n_send*2)) = tgt_coords_z (indices(i))
               buf (ip+(n_send*3)) = dist_dble (i, nb_extra)
            endif
         end do
#endif
!
#ifdef TODO
! TODO: If dest_comp == rank, direct call of psmile_search_donor_extra()
         if (dest_comp == rank) then
!
            search_global%sender = dest
            search_global%msg_extra = msg_extra (1:nd_msg)
!     search_global%len_msg = nd_msg
!
            call psmile_search_donor_extra (search_global, tol, ierror)
            if (ierror > 0) return

         endif
#endif /* TODO */
!
!===> ... Send request to destination process
!         extra_msg is received by routine PSMILe_Enddef_extra and
!         the other buffer are received by routine
!         PSMILe_Enddef_action_extra.
!
!  extra_msg (1)  = Type of request
!  extra_msg (2)  = Datatype
!  extra_msg (3)  = Length of buffer "ibuf" containing the integer data.
!  extra_msg (4)  = Length of buffer  "buf" containing the
!                   real/double data.
!  extra_msg (5)  = Global comp id of component
!  extra_msg (6)  = Global var id "id_trans_out"
!  extra_msg (7)  = Grid Type of sending grid
!                   (PRISM_Reglonlatvrt, ...)
!  extra_msg (8)  = Number of control volumes sent
!  extra_msg (9)  = Number of integer data items per control
!                   volume sent
!  extra_msg (10) = Number of real/double data items per control
!                   volume sent
!  extra_msg (11) = Is data on global offset available ?
!  extra_msg (12) = I-Index of global offset (partition)
!  extra_msg (13) = J-Index of global offset (partition)
!  extra_msg (14) = K-Index of global offset (partition)
!
!  extra_msg (15) = Number of the request
!  extra_msg (16) = num_neigh (Number of interpolation bases)
!  extra_msg (17) = local grid id on message destination process
!
         extra_msg (1)  = PSMILe_nnghbr3D
         extra_msg (2)  = PRISM_DOUBLE_PRECISION
         extra_msg (3)  = len_ibuf
         extra_msg (4)  = len_rbuf
         extra_msg (5)  = comp_info%global_comp_id
         extra_msg (6)  = search%msg_intersections%transient_out_id
!
         extra_msg (7)  = grid_type
         extra_msg (8)  = n_send
         extra_msg (9)  = len_item
         extra_msg (10) = len_rtem
!
!        ... No partition info
!        Partition info senden; dann koennte man srcloc ausnutzen
!
         extra_msg (11) = 0
!
         extra_msg (15) = n_recv_req
         extra_msg (16) = nb_extra

         extra_msg (17) = comp_info%all_extent_infos(1,igrid)
!
#ifdef DEBUGX
      print *, 'extra_msg (1:4)', extra_msg (1:4)
         do i = 1, n_send
            print *, ibuf ((i-1)*len_item+1:i*len_item)
         end do
#endif
!
         call psmile_bsend (extra_msg, nd_msgextra, MPI_INTEGER, &
                            dest, exttag, comm_psmile, ierror)
         if (ierror /= MPI_SUCCESS) then
            ierrp (1) = ierror
            ierrp (2) = dest
            ierrp (3) = exttag

            ierror = PRISM_Error_Send

            call psmile_error (ierror, 'psmile_bsend(msg)', &
                               ierrp, 3, __FILE__, __LINE__ )
            return
         endif
!
         call psmile_bsend (ibuf, len_ibuf, MPI_INTEGER, &
                            dest, exttag, comm_psmile, ierror)
         if (ierror /= MPI_SUCCESS) then
            ierrp (1) = ierror
            ierrp (2) = dest
            ierrp (3) = exttag

            ierror = PRISM_Error_Send

            call psmile_error (ierror, 'psmile_bsend(ibuf)', &
                               ierrp, 3, __FILE__, __LINE__ )
            return
         endif
!
         call psmile_bsend (buf, len_rbuf, MPI_DOUBLE_PRECISION, &
                            dest, exttag, comm_psmile, ierror)
         if (ierror /= MPI_SUCCESS) then
            ierrp (1) = ierror
            ierrp (2) = dest
            ierrp (3) = exttag

            ierror = PRISM_Error_Send

            call psmile_error (ierror, 'psmile_bsend(buf)', &
                               ierrp, 3, __FILE__, __LINE__ )
            return
         endif
!
!===> ... Setup request for answer of destination process
!         ??? wirklich; sollte wohl von dem normalen abgehandelt werden
!
         if (n_recv_req == 1) then
            call MPI_Irecv (answer, nd_msgextra, MPI_INTEGER, MPI_ANY_SOURCE, &
                            rexttag, comm_psmile, recv_req, ierror)
            if (ierror /= MPI_SUCCESS) then

               ierrp (1) = ierror
               ierrp (2) = dest
               ierrp (3) = rexttag

               ierror = PRISM_Error_Recv

               call psmile_error ( ierror, 'MPI_Irecv', &
                                   ierrp, 3, __FILE__, __LINE__ )
               return

            endif
         endif
!
      end do ! igrid
!
!===> Deallocate
!
      Deallocate (send_mask, boxes, igrid_to_dest_comp)
!
      if (ndrbuf > 0) Deallocate (buf)
      if (ndibuf > 0) Deallocate (ibuf)
!
#ifdef __BUG__
      if (n_recv_req == 0) then

!rr      send_info%nrecv and send_info%num2recv already may contain the results
!rr      from the global search for filling the interpolation stencil. Here we
!rr      perform the extra search for points for which no unmasked source points
!rr      have been found. It must not be set to zero here as it is later used
!rr      in psmile_compact_neighbours to compact points from the global search.
!rr      send_info%nrecv and send_info%num2recv get updated in psmile_select_nn_found
!rr      which is called from this routine.

         send_info%nrecv    = 0
         send_info%num2recv = 0

         go to 1000
      endif
#endif

      if (n_recv_req == 0) go to 1000

!
!-------------------------------------------------------------------------
!     Wait for answers
!      Note: (2) Message from a partner with an intersection is disabled
!                since request for grid receive (lrequest (3)) is setup
!                if such a message is found.
!            ??? Koennte man enablen wenn lrequest(3) == MPI_REQUEST_NULL
!            (3) Old receive of grid data (paction%lrequest (3)) is disabled
!-------------------------------------------------------------------------
!
! nrecv      = Number of significant messages received
!
      nrecv = 0
!
      Allocate (sel_info (n_recv_req), stat = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = n_recv_req

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'sel_info', ierrp, 2, &
                             __FILE__, __LINE__ )
         return
      endif
!
!===> Save old values of requests (2:3),
!     ignore index 2 () and
!     set index 3 to first receive of returned meta data
!
      save_lreq (2:3) = paction%lrequest (2:3)
      paction%lrequest (2) = MPI_REQUEST_NULL
!
      do n = 1, n_recv_req
         paction%lrequest (3) = recv_req
!
         index = 0
         do while (index /= 3) 
#ifdef DEBUG
            print *, trim(ch_id), paction%nreq, recv_req
            call psmile_flushstd
#endif
!
            call MPI_Waitany (paction%nreq, paction%lrequest, index, status, ierror)

            if ( ierror /= MPI_SUCCESS ) then
               ierrp (1) = ierror
               ierrp (2) = status (MPI_SOURCE)
               ierrp (3) = status (MPI_TAG)

               ierror = PRISM_Error_MPI

               call psmile_error ( ierror, 'MPI_Waitany', &
                                 ierrp, 3, __FILE__, __LINE__ )
               return
            endif

#ifdef PRISM_ASSERTION
            if (index == MPI_UNDEFINED) then
               call psmile_assert ( __FILE__, __LINE__, &
                                    'request list is empty')
            endif
#endif
!
            if (index /= 3) then
               call psmile_enddef_action (search, index, status, ierror)
               if (ierror > 0) return
            endif
         end do ! while
!
!===> ... Get 
!
         sender = status (MPI_SOURCE)
         len_ibuf = answer (3)
         len_rbuf = answer (4)
!
#ifdef VERBOSE
         print 9960, trim(ch_id), sender, len_ibuf, len_rbuf

         call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
         if (len_ibuf < 0) then
            print *, 'len_ibuf =', len_ibuf
            call psmile_assert (__FILE__, __LINE__, &
                 "len_ibuf should be >= 0")
         endif
!
         if (len_rbuf < 0) then
            print *, 'len_rbuf =', len_rbuf
            call psmile_assert (__FILE__, __LINE__, &
                 "len_rbuf should be >= 0")
         endif
#endif
!
         if (len_ibuf > 0) then
            Allocate (ibuf (1:len_ibuf), stat = ierror)

            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = len_ibuf

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'ibuf', ierrp, 2, &
                                   __FILE__, __LINE__ )
               return
            endif
!
            call MPI_Recv (ibuf, len_ibuf, MPI_INTEGER, sender, &
                           rexttag, comm_psmile, status, ierror)
            if (ierror /= MPI_SUCCESS) then

               ierrp (1) = ierror
               ierrp (2) = sender
               ierrp (3) = rexttag

               ierror = PRISM_Error_Recv

               call psmile_error ( ierror, 'MPI_Recv(ibuf)', ierrp, 3, &
                                   __FILE__, __LINE__ )
               return
            endif
!
            Allocate (buf (1:len_rbuf), stat = ierror)

            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = len_rbuf

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'buf', ierrp, 2, &
                                   __FILE__, __LINE__ )
               return
            endif
!
            call MPI_Recv (buf, len_rbuf, MPI_DOUBLE_PRECISION, sender, &
                           rexttag, comm_psmile, status, ierror)
            if (ierror /= MPI_SUCCESS) then

               ierrp (1) = ierror
               ierrp (2) = sender
               ierrp (3) = rexttag

               ierror = PRISM_Error_Recv

               call psmile_error ( ierror, 'MPI_Recv(rbuf)', ierrp, 3, &
                                   __FILE__, __LINE__ )
               return
            endif
!
!===> ...
!
! n_send  = Number of cells to be searched
! n_found = Number of interpolation bases found
! n_liste = Number of data items to be used for interpolation
!           n_liste <= n_found since a data item can be used
!           multiply for interpolation.
!
            n_send  = answer (7)
            n_liste = answer (8)
            n_found = answer (9)
            ireq    = answer (15)
!
#ifdef PRISM_ASSERTION
            if (answer (1) /= PSMILe_nnghbr3D) then
               print *, 'answer(1)', answer(1), PSMILe_nnghbr3D
               call psmile_assert (__FILE__, __LINE__, &
                    "expected nearest neighbour interpolation")
            endif
!
            if (ireq < 1 .or. ireq > n_recv_req) then
               print *, 'ireq, n_recv_req', ireq, n_recv_req
               call psmile_assert (__FILE__, __LINE__, &
                    "ireq must be in range of 1:n_recv_req")
            endif
#endif
!
            nrecv = nrecv + 1
!
            sel_info(nrecv)%sender    = sender
            sel_info(nrecv)%n_liste   = n_liste
            sel_info(nrecv)%index     = answer (10)
            sel_info(nrecv)%method_id = answer (11)
            sel_info(nrecv)%msg_id    = answer (16)

            sel_info(nrecv)%dble_buf => buf
!
!===> ... Allocate and initialize array "selected"
!
! selected = Temporary vector for global nearest neighbour search
!            for extra points with
!            selected (1, *) = Number (index) of the receive
!                              from which the nearest neigbour was
!                              selected.
!            selected (2, *) = Index in list of points 
!                              (sent by neighbouring process)
!            (cf. routine PSMILe_Add_nn_found_dble)
!            Dimension : selected (2, n_extra)
!
!
            if (.not. Allocated (selected)) then
               Allocate (selected(2, n_extra), stat = ierror)
!
               if (ierror /= 0) then
                  ierrp (1) = n_extra * 2

                  ierror = PRISM_Error_Alloc

                  call psmile_error ( ierror, 'selected', ierrp, 1, &
                              __FILE__, __LINE__ )
                  return
               endif
!
               selected (1, :) = 0
            endif
!
!===> ... Add information on points found by global search 
!
            ip_dist = n_liste * ndim_3d
!
            call psmile_add_nn_found_dble (search, extra_search,     &
                        ibuf (1:n_send),                             &
                        ibuf   (n_send+1:2*n_send), n_send,          &
                        ibuf (2*n_send+1:2*n_send+n_found),          &
                        buf (ip_dist+1:ip_dist+n_found), n_found,    &
                        nb_extra, selected, sel_info, nrecv, ierror)
            if (ierror > 0) return
!
            Deallocate (ibuf)
         endif
!
!===> ... Setup new request
!
         if (n < n_recv_req) then
            call MPI_Irecv (answer, nd_msgextra, MPI_INTEGER, MPI_ANY_SOURCE, &
                              rexttag, comm_psmile, recv_req, &
                              ierror)
            if (ierror /= MPI_SUCCESS) then

               ierrp (1) = ierror
               ierrp (2) = dest
               ierrp (3) = rexttag

               ierror = PRISM_Error_Recv

               call psmile_error ( ierror, 'MPI_Irecv', &
                                    ierrp, 3, __FILE__, __LINE__ )
               return
            endif
         endif
      end do
!
!===> Determine which points are selected and return info to
!     neighbouring processes
!
      if (nrecv > 0) then
         call psmile_select_nn_found (search, extra_search,  &
                     send_info,                              &
                     selected, sel_info, nrecv, nb_extra,    &
                     neighbors_3d, nloc, num_neigh,          &
                     Grids(grid_id)%grid_shape, ierror)
         if (ierror > 0) return
      endif
!
!    ... Free memory
!
      if (Allocated (selected)) then
         Deallocate (selected)
      endif
!
      Deallocate (sel_info)
!
!    ... Restore original requests (2:3)
!
#ifdef PRISM_ASSERTION
      if (paction%lrequest (2) /= MPI_REQUEST_NULL .or. &
          paction%lrequest (3) /= MPI_REQUEST_NULL) then
         print *, 'request: ', paction%lrequest (2:3)
         call psmile_assert ( __FILE__, __LINE__, &
                             'Illegal request stored')

      endif
#endif
!
      paction%lrequest (2:3) = save_lreq (2:3)
!
!-------------------------------------------------------------------------
!===> All done
!-------------------------------------------------------------------------
!
1000  continue
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

#ifdef VERBOSE

9990 format (1x, a, ': psmile_global_search_nnx_dble: var_id', i3, &
                    ' to ', i3, '(', i2, ')')
9980 format (1x, a, ': psmile_global_search_nnx_dble: eof ierror =', i3)
9970 format (1x, a, ': psmile_global_search_nnx_dble: send from', i3, &
                    ' to', i3, '[', i3, '], n_send =', i6)

9960 format (1x, a, ': psmile_global_search_nnx_dble: got rexttag message:', &
                    ' sender ', i4, ', len_ibuf, len_rbuf', 2i8)
9950 format (1x, a, ': psmile_global_search_nnx_dble: before waitany :', &
                    'nreq =', i4, ', recv_req ', i8)
#endif /* VERBOSE */

#ifdef DEBUG
#endif

      end subroutine PSMILe_global_search_nnx_dble