psmile_global_search_cell_real.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! $COPYRIGHT_NLE-IT
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_global_search_cell_real
!
! !INTERFACE:

      subroutine psmile_global_search_cell_real ( grid_id, var_id,  &
                  comp_info, send_info, search, extra_search, ncpl, &
                  nbr_cells, n_intmethods, interpolation_methods,   &
                  interpolation_search, ierror )
!
! !USES:
!
      use PRISM
!
      use PSMILe, dummy_interface => PSMILe_global_search_cell_real

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent(In)                  :: grid_id

!     grid id of source grid

      Integer, Intent(In)                  :: var_id

!     field id of source data

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

      Type (Send_information), Intent (InOut) :: send_info
!
!     Metadata for target locations
!
      Type (Enddef_search), Intent (InOut) :: search
!
!     Info's on coordinates to be searched
!
      Type (Extra_search_info), Intent (InOut) :: extra_search
!
!     Info's on global coordinates to be searched
!
      Integer, Intent (In)                 :: ncpl
!
!     Number of target cells treated on this source process
!
      Integer, Intent (InOut)              :: nbr_cells(ncpl)
!
!     Number of source cells per target cell
!
      Integer, Intent (In)                 :: n_intmethods
!
!     Number of interpolations methods used to interpolate the data.
!     Dimension of vectors "interpolation_methods" and "interpolation_search".

      Integer, Intent (In)                 :: interpolation_methods (n_intmethods)
!
!     Codes for the interpolation methods
!
      Logical, Intent (In)                 :: interpolation_search  (n_intmethods)
!
!     Is global search required for corresponding interpolation method ?
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)                :: ierror

!     Returns the error code of PSMILe_global_search_cell_real;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS
!
!     ... for send to partners
!
      Integer, Parameter                 :: n_extra = 7
      Integer, Parameter                 :: send_appl_index = 5

      Logical, Pointer                   :: lbuf (:)
      Real, Pointer                      :: buf (:)

      Type (real_vector), Pointer        :: real_buf(:)
      Type (integer_vector), Allocatable :: tgt_ibuf(:)
      Type (integer_vector), Allocatable :: src_ibuf(:)
      Type (logical_vector), Allocatable :: log_buf(:)

      Type (Corner_Block), Pointer       :: cp
      Type (Method),       Pointer       :: mp

      Integer, Pointer                   :: new_neighcells_3d(:,:,:)
      Integer, Pointer                   :: neigh_3d(:,:)
      Logical, Pointer                   :: mask_array(:)
      Integer                            :: mask_shape (2, ndim_3d)

      Integer, Allocatable               :: rextra_msg (:,:)
      Integer, Allocatable               :: nbr_tgt_cells(:)
      Integer, Allocatable               :: nbr_src_cells(:)
      Integer, Allocatable               :: rank_list(:)
      Integer, Allocatable               :: cell_trans(:,:)
      Integer, Allocatable               :: save_lreq(:)
      Integer, Allocatable               :: lrequest(:)
      Integer, Allocatable               :: tgt_list(:), tgt_cell(:), tgt_index(:)

      Logical                            :: src_mask_available

      Integer                            :: tgt_idx
      Integer                            :: new_nbr_cells(ncpl)
      Integer                            :: mask_id
      Integer                            :: method_id
      Integer                            :: nbr_corners
      Type (enddef_msg_locations)        :: msg_locations
      Integer                            :: msgloc (msgloc_size)
      Integer                            :: extra_msg (n_extra)
      Integer                            :: len
      Integer                            :: lenx, leny, lenxy, m_lenx
      Integer                            :: rank, idx1, idx2, idx3, idx4, stride
      Integer                            :: nbr_procs(2)
      Integer                            :: interp
      Integer                            :: dest
      Integer                            :: sender
      Integer                            :: status(MPI_STATUS_SIZE)
      Integer                            :: nbr_cells_tot
      Integer                            :: nreq
      Integer                            :: recv_req
      Integer                            :: nrecv, nsend
      Integer                            :: index
      Integer                            :: nc
      Integer                            :: l
      Integer                            :: m, mm
      Integer                            :: n, nn
!
!     ... for error handling
!
      Integer, Parameter                 :: nerrp = 3
      Integer                            :: ierrp (nerrp)
!
!----------------------------------------------------------------------
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_global_search_cell_real" collects source cells
! for target cells located on process boundaries of the source grid.
! Corresponding MPI calls are located in psmile_enddef_action_cell.
!
!   - exchange boundary cell info by sending back
!     information to the source grid processes holding the
!     respective target cells.
!   - determine who is master source for the target cells
!   - let the masters collect the data
!   - store extra source cells on the master source processes
!   - delete real entries on the non-master sides
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 13.08.08    R. Redler    created
! 08.11.06    R. Redler    revised
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_global_search_cell_real.F90 2887 2011-01-14 09:25:49Z redler $
!  $Author: redler $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_global_search_cell_real.F90 2887 2011-01-14 09:25:49Z redler $'
!
!----------------------------------------------------------------------
!  1.) Initialization
!----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9990, trim(ch_id), grid_id, search%msg_intersections%first_tgt_method_id, search%sender

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0

      method_id = Fields(var_id)%method_id
      mask_id = Fields(var_id)%mask_id

      cp => Grids(grid_id)%corner_pointer
      mp => Methods(method_id)
      nreq = paction%nreq

      src_mask_available = mask_id /= PRISM_UNDEFINED

      if (src_mask_available) then
#ifdef DEBUG
        print *, ' Source mask is available.'
#endif  /* DEBUG */
         mask_array => Masks(mask_id)%mask_array
         mask_shape =  Masks(mask_id)%mask_shape
      endif

      new_nbr_cells = nbr_cells
      send_info%nloc = ncpl
      nbr_cells_tot  = sum(nbr_cells)

      if ( Grids(grid_id)%grid_type == PRISM_Irrlonlat_regvrt ) then

         nbr_corners = cp%nbr_corners/2

      else if ( Grids(grid_id)%grid_type == PRISM_Reglonlatvrt ) then

         nbr_corners = 2

      else if ( Grids(grid_id)%grid_type == PRISM_Gaussreduced_regvrt ) then

         nbr_corners = 2

      else

         ierrp (1) = search%grid_type
         ierror = PRISM_Error_Internal
         call psmile_error ( ierror, 'unsupported search grid type', &
              ierrp, 1, __FILE__, __LINE__ )
         return

      endif

      Allocate(save_lreq(nreq), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = nreq
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'save_lreq', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

#ifdef PRISM_ASSERTION
      if (interpolation_methods(1) /= PSMILe_conserv2D ) then
         call psmile_assert (__FILE__, __LINE__, &
              "No global search for this interpolation method")
      endif
      if ( n_intmethods > 1 ) then
         call psmile_assert (__FILE__, __LINE__, &
              "2nd (vertical) interpolation method not yet supported.")
      endif
#endif /* PRISM_ASSERTION */
      if (count(interpolation_search(:)) > 1) then
         ierrp (1) = count(interpolation_search(:))

         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'Global search supported only for one interpolation method', &
              ierrp, 1, __FILE__, __LINE__ )

         return
      endif

      do interp = 1, n_intmethods
         if (interpolation_search(interp) ) exit
      end do

      if (interp > n_intmethods) then
         ierror = PRISM_Error_Internal
         ierrp (1) = n_intmethods

         call psmile_error ( ierror, 'No Global search for interpolation methods', &
              ierrp, 1, __FILE__, __LINE__ )
         return
      endif

      len = 0

      if ( associated(search%boundary_cell) ) then

         len = size(search%boundary_cell(:,1))

         Allocate(cell_trans(len,3), stat=ierror)

         if (ierror /= 0) then
            ierrp (1) = 3*len
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'cell_trans', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

      endif
!
!-------------------------------------------------------------------------
! 2.) Exchange data between source processes
!-------------------------------------------------------------------------
!
!     Send back information about target cells located on source boundary
!     to target process
!
!===> ... Send request to destination process
!         msgloc is received by routine
!         PSMILe_Enddef_action.
!
!
      call psmile_init_enddef_msg_locs (msg_locations)

      msg_locations%requires_conserv_remap = interpolation_methods (interp)
      msg_locations%src_rank               = psmile_rank
      msg_locations%msg_len                = len
      msg_locations%transi_in_id           = search%msg_intersections%transient_in_id
      msg_locations%transi_out_id          = search%msg_intersections%transient_out_id

      call psmile_pack_msg_locations (msg_locations, msgloc)

      dest = search%sender

#ifdef DEBUG
      print *, ' Bsend tag ', loctag+search%msg_intersections%relative_msg_tag, ' to ', dest
#endif /* DEBUG */
      call psmile_bsend ( msgloc, msgloc_size, MPI_INTEGER, dest,           &
                          loctag+search%msg_intersections%relative_msg_tag, &
                          comm_psmile, ierror )

      if ( ierror /= MPI_SUCCESS ) then
         ierrp (1) = ierror
         ierrp (2) = dest
         ierrp (3) = loctag+search%msg_intersections%relative_msg_tag
         ierror = PRISM_Error_Send
         call psmile_error ( ierror, 'psmile_bsend(msg)', &
                             ierrp, 3, __FILE__, __LINE__ )
         return
      endif
      !
      ! Transfer indices and information about those target cells
      !  which touch the boundary of this source domain, received
      !  by routine PSMILe_Enddef_action_cell.
      !
      if ( len > 0 ) &
      call psmile_bsend ( search%boundary_cell, 5*len, MPI_INTEGER, &
                          dest, celltag, comm_psmile, ierror )

      if ( ierror /= MPI_SUCCESS ) then
         ierrp (1) = ierror
         ierrp (2) = dest
         ierrp (3) = celltag
         ierror = PRISM_Error_Send
         call psmile_error ( ierror, 'psmile_bsend(msg)', &
                             ierrp, 3, __FILE__, __LINE__ )
         return
      endif
      !
      ! Receive the number of source processes to which the source
      ! has to transfer data, sent by routine PSMILe_Enddef_action_cell.
      !
      call MPI_Irecv ( nbr_procs, 2, MPI_INTEGER, &
                       dest, celltag+1, comm_psmile, recv_req, ierror )

      if ( ierror /= MPI_SUCCESS ) then
         ierrp (1) = ierror
         ierrp (2) = dest
         ierrp (3) = celltag+1
         ierror = PRISM_Error_Recv
         call psmile_error ( ierror, 'MPI_Recv', &
                             ierrp, 3, __FILE__, __LINE__ )
         return
      endif

#ifdef DEBUG
      print *, ' Posted Irecv request ', recv_req, ' from ', dest
#endif /* DEBUG */

      save_lreq = paction%lrequest
      !
      ! Give absolute priority to cell requests by temporarily deactivating
      ! all other pending request.
      !
      ! TODO: For the sake of performance it has to be checked whether there
      !       are any requests which can savely be handled in between without
      !       running into a deadlock.
      !
      paction%lrequest = MPI_REQUEST_NULL

      paction%lrequest (3) = recv_req

!     Are there any requests that we should better deactivate temporaily?
!
!     paction%lrequest (?:?) = MPI_REQUEST_NULL
!
      index = 0

      do while ( index /= 3)

#ifdef DEBUG
         do n = 1, nreq
           print *, ' paction%lrequest ', n, paction%lrequest(n)
         enddo
#endif /* DEBUG */

         call MPI_Waitany ( nreq, paction%lrequest, index, status, ierror )

#ifdef DEBUG
         print *, ' Processing with index ', index, ' from ', status(MPI_SOURCE)
#endif /* DEBUG */

         if ( index /= 3 ) then

            call psmile_enddef_action (search, index, status, ierror )

         endif

      enddo

      paction%lrequest(3) = MPI_REQUEST_NULL

      nsend = nbr_procs(1)
      nrecv = nbr_procs(2)

#ifdef DEBUG
      print *, ' This rank has to serve        ', nsend, ' source grid processes.'
      print *, ' This rank has to receive from ', nrecv, ' source grid processes.'
#endif /* DEBUG */

      if ( max(nsend,nrecv) > 0 ) then
         Allocate ( tgt_ibuf (max(nsend,nrecv)), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = max(nsend,nrecv)
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'tgt_ibuf', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif
      endif

      if ( nsend > 0 ) then
         Allocate(real_buf      (nsend), &
                  log_buf       (nsend), &
                  src_ibuf      (nsend), &
                  nbr_tgt_cells (nsend), &
                  nbr_src_cells (nsend), &
                  rank_list     (nsend), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 6*nsend
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'nsend buffer', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         nbr_tgt_cells = 0
         nbr_src_cells = 0
         rank_list = -1
      endif

      !
      ! Receive the list of target indices and related information from
      ! target routine PSMILe_Enddef_action_cell for which source cells
      ! have to be transferred to other source processes.
      ! cell_trans target cell index on the sender source contains a
      ! sorted list with increasing order ranks stored in
      ! cell_trans(:,2)
      !
      ! cell_trans(:,1) : target cell index on the sender source
      ! cell_trans(:,2) : rank of receiver source process
      ! cell_trans(:,3) : target cell index on the receiver source
      !

#ifdef DEBUGX
      print *, ' Receiving data ', len, ' records from ', dest
      if ( len == 0 ) &
      print *, ' Receive canceled because of zero message '
#endif /* DEBUGX */
      if ( len > 0 ) &
      call MPI_Recv ( cell_trans, 3*len, MPI_INTEGER, &
                      dest, celltag+2, comm_psmile, status, ierror )

      if ( ierror /= MPI_SUCCESS ) then
         ierrp (1) = ierror
         ierrp (2) = dest
         ierrp (3) = celltag+2
         ierror = PRISM_Error_Recv
         call psmile_error ( ierror, 'MPI_Recv', &
                             ierrp, 3, __FILE__, __LINE__ )
         return
      endif
      !
      ! Determine number of data to be transferred
      !   to each destination source
      !
      if ( nsend > 0 ) then

#ifdef DEBUGX
         do n = 1, len
!            if ( cell_trans(n,1) > -1 ) then
               print 9970, trim(ch_id), n, cell_trans(n,1), &
                    ' -> ', cell_trans(n,2), cell_trans(n,3)
!            endif
         enddo
9970     format (1x, a, ': ', i6, ' Transfer index ', i6, a4, i3, i6)
#endif /* DEBUGX */

         nn = 1
         m  = 0

         rank_list(nn) = cell_trans(1,2)

         do n = 1, len
            if ( n > 1 ) then
               if ( cell_trans(n,2) > cell_trans(n-1,2) ) then
                  nn = nn + 1
                  rank_list(nn) = cell_trans(n,2)
               endif
            endif
            if ( cell_trans(n,2) == -1 ) exit

            tgt_idx = cell_trans(n,1) ! search%boundary_cell(n,1) ! local index of target cell
            nbr_tgt_cells(nn) = nbr_tgt_cells(nn) + 1
            nbr_src_cells(nn) = nbr_src_cells(nn) + nbr_cells(tgt_idx)
         enddo
#ifdef DEBUG
         do n = 1, nsend
            print *, nbr_tgt_cells(n), ' target cell info to transfer.'
            print *, nbr_src_cells(n), ' source cells have to be transferred to rank ', rank_list(n)
         enddo
#endif /* DEBUG */

!
!-------------------------------------------------------------------------
! 3.) Send data to the designated child source processes
!-------------------------------------------------------------------------
!
! TODO: Eliminate multiple entries from cell_trans array and create extra index.
!
!
         do n = 1, nsend
            Allocate ( tgt_ibuf(n)%vector(nbr_tgt_cells(n)*2), &
                       src_ibuf(n)%vector(nbr_src_cells(n)*3), &
                       log_buf(n)%vector(nbr_src_cells(n)),    &
                       real_buf(n)%vector(nbr_src_cells(n)*2*nbr_corners), stat=ierror )
            if (ierror /= 0) then
               ierrp (1) = 2*nbr_tgt_cells(n) + (2*nbr_corners+4)*nbr_src_cells(n)
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'src/tgt_buf(n)%ibuf/real_buf ', &
                  ierrp, 1, __FILE__, __LINE__ )
               return
            endif
         enddo
         !
         ! Fill in target index information
         !
         nn = 1
         m  = 0

         do n = 1, len

            if ( n > 1 ) then
               if ( cell_trans(n,2) > cell_trans(n-1,2) ) then
                  nn = nn + 1
                  m = 0
               endif
            endif

            if ( cell_trans(n,2) == -1 ) exit

            tgt_idx = cell_trans(n,1) ! search%boundary_cell(n,1) ! local index of target cell
            m = m+1

            tgt_ibuf(nn)%vector(2*m-1) = cell_trans(n,3)
            tgt_ibuf(nn)%vector(2*m  ) = nbr_cells(tgt_idx)

         enddo

         !
         ! Distribute information to the different send buffers
         !

         select case (Grids(grid_id)%grid_type)

         case (PRISM_Irrlonlat_regvrt)

            !
            ! Irrlonlat_regvrt source grids
            !

            lenx  = (cp%corner_shape(2,1)-cp%corner_shape(1,1)+1)
            lenxy = (cp%corner_shape(2,2)-cp%corner_shape(1,2)+1) * lenx

            if (src_mask_available) then
               m_lenx = (mask_shape(2,1)-mask_shape(1,1)+1)
#ifdef PRISM_ASSERTION
               if ( cp%corner_shape(1,2) /= mask_shape(1,2) .or. &
                    cp%corner_shape(1,1) /= mask_shape(1,1) ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "Mask and corners have to be of the same shape.")
               endif
               if ( mask_shape(2,3) - mask_shape(1,3) + 1 > 1 ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "3d masks not yet supported.")
               endif
#endif /* PRISM_ASSERTION */
            endif

            l  = 0
            rank = 1

            do n = 1, len

               if ( n > 1 ) then
                  if ( cell_trans(n,2) > cell_trans(n-1,2) ) then
                     rank = rank + 1
                     l    = 0
                  endif
               endif
               if ( cell_trans(n,2) == -1 ) exit

               stride = nbr_corners * nbr_src_cells(rank)

               tgt_idx = cell_trans(n,1) ! search%boundary_cell(n,1)
#ifdef DEBUGX
               print *, ' Corners for local target index ', tgt_idx
#endif /* DEBUGX */
               nn = sum(nbr_cells(1:tgt_idx-1))

               do mm = 1, nbr_cells(tgt_idx)

                  idx3 = (neighcells_3d(2,nn+mm,1)-cp%corner_shape(1,2))*lenx + &
                          neighcells_3d(1,nn+mm,1)-cp%corner_shape(1,1)+1

                  src_ibuf(rank)%vector(3*(l+mm)-2) = neighcells_3d(1,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)-1) = neighcells_3d(2,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)  ) = neighcells_3d(3,nn+mm,1)

                  do nc = 1, nbr_corners
                     idx1 = (nc-1)*nbr_src_cells(rank)+l+mm
                     idx2 = stride + idx1
                     real_buf(rank)%vector(idx1) = cp%corners_real(1)%vector((nc-1)*lenxy+idx3)
                     real_buf(rank)%vector(idx2) = cp%corners_real(2)%vector((nc-1)*lenxy+idx3)
#ifdef DEBUGX
                     print *, ' send corners ', l+mm, real_buf(rank)%vector(idx1), &
                                                      real_buf(rank)%vector(idx2)
                  end do

                  print *, ' 2d indices   ', neighcells_3d(1,nn+mm,1), &
                                             neighcells_3d(2,nn+mm,1)
#else
                  end do
#endif /* DEBUGX */
               enddo ! mm

               if (src_mask_available) then

                  do mm = 1, nbr_cells(tgt_idx)
                     idx3 = (neighcells_3d(2,nn+mm,1)-mask_shape(1,2))*m_lenx + &
                             neighcells_3d(1,nn+mm,1)-mask_shape(1,1)+1
                     log_buf(rank)%vector(l+mm) = mask_array(idx3)
                  enddo
#ifdef DEBUGX
                  print *, ' send mask is ', log_buf(rank)%vector(l+1:l+nbr_cells(tgt_idx))
#endif
               endif

               ! Erase target cell entries and related info which
               ! have been transfered to a remote source process
#ifdef DEBUGX
               print *, ' Erase target cell index ', tgt_idx
               print *
#endif /* DEBUGX */
               new_nbr_cells(tgt_idx) = 0
               send_info%nloc = send_info%nloc - 1

               ! Increment indices

               l = l + nbr_cells(tgt_idx)

            enddo ! len

         case (PRISM_Reglonlatvrt)

            !
            ! Reglonlatvrt source grids
            !

            lenx = (cp%corner_shape(2,1)-cp%corner_shape(1,1)+1)
            leny = (cp%corner_shape(2,2)-cp%corner_shape(1,2)+1)

            if (src_mask_available) then
               m_lenx = (mask_shape(2,1)-mask_shape(1,1)+1)
#ifdef PRISM_ASSERTION
               if ( cp%corner_shape(1,2) /= mask_shape(1,2) .or. &
                    cp%corner_shape(1,1) /= mask_shape(1,1) ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "Mask and corners have to be of the same shape.")
               endif
               if ( mask_shape(2,3) - mask_shape(1,3) + 1 > 1 ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "3d masks not yet supported.")
               endif
#endif /* PRISM_ASSERTION */
            endif

            l  = 0
            rank = 1

            do n = 1, len

               if ( n > 1 ) then
                  if ( cell_trans(n,2) > cell_trans(n-1,2) ) then
                     rank = rank + 1
                     l    = 0
                  endif
               endif
               if ( cell_trans(n,2) == -1 ) exit

               stride = nbr_corners * nbr_src_cells(rank)

               tgt_idx = cell_trans(n,1) ! search%boundary_cell(n,1)
#ifdef DEBUGX
               print *, ' Corners for local target index ', tgt_idx
#endif /* DEBUGX */
               nn = sum(nbr_cells(1:tgt_idx-1))

               do mm = 1, nbr_cells(tgt_idx)

                  src_ibuf(rank)%vector(3*(l+mm)-2) = neighcells_3d(1,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)-1) = neighcells_3d(2,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)  ) = neighcells_3d(3,nn+mm,1)

                  idx3 = neighcells_3d(1,nn+mm,1)-cp%corner_shape(1,1)+1
                  idx4 = neighcells_3d(2,nn+mm,1)-cp%corner_shape(1,2)+1

                  do nc = 1, nbr_corners
                     idx1 = (nc-1)*nbr_src_cells(rank)+l+mm
                     idx2 = stride + idx1
                     real_buf(rank)%vector(idx1) = cp%corners_real(1)%vector((nc-1)*lenx+idx3)
                     real_buf(rank)%vector(idx2) = cp%corners_real(2)%vector((nc-1)*leny+idx4)
#ifdef DEBUGX
                     print *, ' send corners ', l+mm, real_buf(rank)%vector(idx1), &
                                                      real_buf(rank)%vector(idx2)
                  end do

                  print *, ' 2d indices   ', neighcells_3d(1,nn+mm,1), &
                                             neighcells_3d(2,nn+mm,1)
#else
                  end do
#endif /* DEBUGX */
               enddo ! mm

               if (src_mask_available) then

                  do mm = 1, nbr_cells(tgt_idx)
                     idx3 = (neighcells_3d(2,nn+mm,1)-mask_shape(1,2))*m_lenx + &
                             neighcells_3d(1,nn+mm,1)-mask_shape(1,1)+1
                     log_buf(rank)%vector(l+mm) = mask_array(idx3)
                  enddo
#ifdef DEBUGX
                  print *, ' send mask is ', log_buf(rank)%vector(l+1:l+nbr_cells(tgt_idx))
#endif
               endif

               ! Erase target cell entries and related info which
               ! have been transfered to a remote source process
#ifdef DEBUGX
               print *, ' Erase target cell index ', tgt_idx
               print *
#endif /* DEBUGX */
               new_nbr_cells(tgt_idx) = 0
               send_info%nloc = send_info%nloc - 1

               ! Increment indices

               l = l + nbr_cells(tgt_idx)

            enddo ! len

         case (PRISM_Gaussreduced_regvrt)

            !
            ! Gauss-reduced source grids
            !

            lenx = (cp%corner_shape(2,1)-cp%corner_shape(1,1)+1)

            if (src_mask_available) then
               m_lenx = (mask_shape(2,1)-mask_shape(1,1)+1)
#ifdef PRISM_ASSERTION
               if ( cp%corner_shape(1,2) /= mask_shape(1,2) .or. &
                    cp%corner_shape(1,1) /= mask_shape(1,1) ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "Mask and corners have to be of the same shape.")
               endif
               if ( mask_shape(2,3) - mask_shape(1,3) + 1 > 1 ) then
                  call psmile_assert (__FILE__, __LINE__, &
                       "3d masks not yet supported.")
               endif
#endif /* PRISM_ASSERTION */
            endif

            l  = 0
            rank = 1

            do n = 1, len

               if ( n > 1 ) then
                  if ( cell_trans(n,2) > cell_trans(n-1,2) ) then
                     rank = rank + 1
                     l    = 0
                  endif
               endif
               if ( cell_trans(n,2) == -1 ) exit

               stride = nbr_corners * nbr_src_cells(rank)

               tgt_idx = cell_trans(n,1) ! search%boundary_cell(n,1)
#ifdef DEBUGX
               print *, ' Corners for local target index ', tgt_idx
#endif /* DEBUGX */
               nn = sum(nbr_cells(1:tgt_idx-1))

               do mm = 1, nbr_cells(tgt_idx)

                  src_ibuf(rank)%vector(3*(l+mm)-2) = neighcells_3d(1,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)-1) = neighcells_3d(2,nn+mm,1)
                  src_ibuf(rank)%vector(3*(l+mm)  ) = neighcells_3d(3,nn+mm,1)

                  idx3 = neighcells_3d(1,nn+mm,1)-cp%corner_shape(1,1)+1
                  idx4 = neighcells_3d(2,nn+mm,1)-cp%corner_shape(1,2)+1

                  do nc = 1, nbr_corners
                     idx1 = (nc-1)*nbr_src_cells(rank)+l+mm
                     idx2 = stride + idx1
                     real_buf(rank)%vector(idx1) = cp%corners_real(1)%vector((nc-1)*lenx+idx3)
                     real_buf(rank)%vector(idx2) = cp%corners_real(2)%vector((nc-1)*lenx+idx4)
#ifdef DEBUGX
                     print *, ' send corners ', l+mm, real_buf(rank)%vector(idx1), &
                                                      real_buf(rank)%vector(idx2)
                  end do

                  print *, ' 2d indices   ', neighcells_3d(1,nn+mm,1), &
                                             neighcells_3d(1,nn+mm,1)
#else
                  end do
#endif /* DEBUGX */
               enddo ! mm

               if (src_mask_available) then

                  do mm = 1, nbr_cells(tgt_idx)
                     idx3 = neighcells_3d(1,nn+mm,1)-mask_shape(1,1)+1
                     log_buf(rank)%vector(l+mm) = mask_array(idx3)
                  enddo
#ifdef DEBUGX
                  print *, ' send mask is ', log_buf(rank)%vector(l+1:l+nbr_cells(tgt_idx))
#endif
               endif

               ! Erase target cell entries and related info which
               ! have been transfered to a remote source process
#ifdef DEBUGX
               print *, ' Erase target cell index ', tgt_idx
               print *
#endif /* DEBUGX */
               new_nbr_cells(tgt_idx) = 0
               send_info%nloc = send_info%nloc - 1

               ! Increment indices

               l = l + nbr_cells(tgt_idx)

            enddo ! len

         case DEFAULT

            ierrp (1) = search%grid_type
            ierror = PRISM_Error_Internal
            call psmile_error ( ierror, 'unsupported search grid type', &
                 ierrp, 1, __FILE__, __LINE__ )
            return

         end select

         nbr_cells_tot = sum(new_nbr_cells)

      endif ! nsend > 0
      !
      ! Reset requests and deallocate memory
      !
!rr   paction%lrequest(2:3) = save_lreq(2:3)

      paction%lrequest = save_lreq

      if ( nsend > nreq ) then
         Deallocate (save_lreq, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = nreq
            ierror = PRISM_Error_Dealloc
            call psmile_error ( ierror, 'save_lreq', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif
      endif
      !
      ! Send data away to the designated parent source processes
      !
      do n = 1, nsend
         !
         ! Prepare header message to be send to destination source process.
         !
         ! (1) = interpolation type
         ! (2) = data type
         ! (3) = number of target cells to be treated
         ! (4) = number of source cell corners
         ! (5) = number of source cells
         !
         extra_msg(1) = interpolation_methods(1)
         extra_msg(2) = PRISM_REAL
         extra_msg(3) = nbr_tgt_cells(n)
         extra_msg(4) = nbr_corners
         extra_msg(5) = nbr_src_cells(n)
         extra_msg(6) = mask_id
         !
         ! This should be sufficient for the later operations during the prism_put phase
         ! -> see psmile_put_field_real and psmile_extract_indices_3d_real
         !
         call psmile_get_info_index (method_id, send_appl_index, idx1, ierror)
         if (ierror > 0) return

         extra_msg(7) = mp%send_infos_appl(idx1)%msg_id

         mp%send_infos_appl(idx1)%dest   = rank_list(n)
         mp%send_infos_appl(idx1)%nloc   = nbr_src_cells(n)

         allocate(neigh_3d(3,nbr_src_cells(n)), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 3*nbr_src_cells(n)
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'neigh_3d', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         do m = 1, nbr_src_cells(n)
            neigh_3d(1,m) = src_ibuf(n)%vector(3*m-2)
            neigh_3d(2,m) = src_ibuf(n)%vector(3*m-1)
            neigh_3d(3,m) = src_ibuf(n)%vector(3*m  )
         enddo

         mp%send_infos_appl(idx1)%neigh_3d => neigh_3d

         call psmile_store_send_info (var_id,                   &
                     search%msg_intersections%transient_out_id, &
                     PRISM_UNDEFINED, PRISM_UNDEFINED, idx1, ierror)
         !
         ! Send header to destination sources
         !
#ifdef DEBUG
         print *, ' Posting send request to psmile rank ', rank_list(n)
#endif /* DEBUG */
         call psmile_bsend ( extra_msg, n_extra, MPI_INTEGER, &
              rank_list(n), celltag+3, comm_psmile, ierror )

        if (ierror /= MPI_SUCCESS) then
            ierrp (1) = ierror
            ierrp (2) = rank_list(n)
            ierrp (3) = celltag+3

            ierror = PRISM_Error_Send

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

#ifdef DEBUGX
         do m = 1, nbr_tgt_cells(n)
            print '(a16,i6,a13,i4)', ' send tgt index ',  tgt_ibuf(n)%vector(2*m-1), &
                                     ' # src cells ',     tgt_ibuf(n)%vector(2*m  )
         enddo
#endif /* DEBUGX */
         call psmile_bsend ( tgt_ibuf(n)%vector, 2*nbr_tgt_cells(n), &
              MPI_INTEGER, rank_list(n), celltag+4, comm_psmile, ierror )

         call psmile_bsend ( real_buf(n)%vector, 2*nbr_corners*nbr_src_cells(n), &
              MPI_REAL, rank_list(n), celltag+6, comm_psmile, ierror )

         if ( src_mask_available ) &
         call psmile_bsend ( log_buf(n)%vector, nbr_src_cells(n), &
              MPI_LOGICAL, rank_list(n), celltag+7, comm_psmile, ierror )
         !
         ! Free send memory
         !
         Deallocate ( tgt_ibuf(n)%vector, src_ibuf(n)%vector, real_buf(n)%vector, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 2*nbr_tgt_cells(n)+(2*nbr_corners+1)*nbr_src_cells(n)
            ierror = PRISM_Error_Dealloc
            call psmile_error ( ierror, 'src/tgt buffer', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

      enddo ! n = 1, nsend

      if ( nsend > 0 ) then
         Deallocate ( real_buf, src_ibuf, nbr_tgt_cells, nbr_src_cells, rank_list, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 5*nsend
            ierror = PRISM_Error_Dealloc
            call psmile_error ( ierror, 'src_ibuf', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif
      endif
!
!-------------------------------------------------------------------------
! 4.) Receive data from the designated child source processes
!-------------------------------------------------------------------------
!
      extra_search%nrecv = 0

      if ( nrecv > 0 ) then

         extra_search%nrecv = nrecv

         Allocate (rextra_msg(n_extra,nrecv),      &
                   extra_search%dble_bufs(nrecv),  &
                   extra_search%log_bufs(nrecv),   &
                   lrequest(nrecv),                &
                   send_info%len_sent(nrecv),      &
                   send_info%sender_global(nrecv), &
                   send_info%msg_id(nrecv),        &
                   stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = (n_extra+4)*nrecv
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'rextra_msg, lrequest, ... ', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         lrequest = MPI_REQUEST_NULL

         do n = 1, nrecv

            call MPI_Irecv ( rextra_msg(:,n), n_extra, MPI_INTEGER, &
                 MPI_ANY_SOURCE, celltag+3, comm_psmile, lrequest(n), ierror )

            if ( ierror /= MPI_SUCCESS ) then
               ierrp (1) = ierror
               ierrp (2) = MPI_ANY_SOURCE
               ierrp (3) = celltag+3

               ierror = PRISM_Error_MPI

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

         enddo

         do n = 1, nrecv

            call MPI_Waitany ( nrecv, 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

            sender = status(MPI_SOURCE)
#ifdef DEBUG
            print *, ' Request index ', index , ' received from rank ', sender
#endif /* DEBUG */

            len = 2*rextra_msg(3,index)

            Allocate ( tgt_ibuf(index)%vector(len), stat = ierror)
            if (ierror /= 0) then
               ierrp (1) = len
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'tgt_ibuf%vector', ierrp, 1, &
                    __FILE__, __LINE__ )
               return
            endif

            ! Receive transferred target cell indices

            call MPI_Recv ( tgt_ibuf(index)%vector, len, MPI_INTEGER, sender, &
                 celltag+4, comm_psmile, status, ierror )
#ifdef DEBUGX
            do m = 1, rextra_msg(3,n)
               print '(a16,i6,a13,i4)', &
                     ' recv tgt index ', tgt_ibuf(index)%vector(2*m-1), &
                     ' # src cells ',    tgt_ibuf(index)%vector(2*m)
            enddo
#endif /* DEBUGX */

            len = 2*rextra_msg(4,index)*rextra_msg(5,index)

            Allocate ( buf(len), stat=ierror)
            if (ierror /= 0) then
               ierrp (1) = len
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'buf', ierrp, 1, &
                    __FILE__, __LINE__ )
               return
            endif

            ! Receive extra source cell coordinates

            call MPI_Recv ( buf, len, MPI_REAL, &
                 sender, celltag+6, comm_psmile, status, ierror )

            extra_search%real_bufs(index)%vector => buf

#ifdef DEBUGX
            stride = rextra_msg(5,index)*rextra_msg(4,index)

            do m = 1, rextra_msg(5,index)
               do nc = 1, rextra_msg(4,index) ! nbr_corners
                  idx1 = (nc-1)*rextra_msg(5,index)+m
                  idx2 = stride + idx1
                  print *, ' receive corners ', m, &
                    extra_search%real_bufs(index)%vector(idx1), &
                    extra_search%real_bufs(index)%vector(idx2)
               end do
               print *
            enddo
#endif /* DEBUGX */

            ! Receive extra mask values

            if ( rextra_msg(6,index) /= PRISM_UNDEFINED ) then

               len = rextra_msg(5,index)
               Allocate ( lbuf(len), stat=ierror)
               if (ierror /= 0) then
                  ierrp (1) = len
                  ierror = PRISM_Error_Alloc
                  call psmile_error ( ierror, 'buf', ierrp, 1, &
                       __FILE__, __LINE__ )
                  return
               endif
#ifdef DEBUGX
               print *, ' receive extra mask values '
#endif
               call MPI_Recv ( lbuf, len, MPI_LOGICAL, &
                      sender, celltag+7, comm_psmile, status, ierror )

#ifdef DEBUGX
               print *, ' received extra mask values '
#endif
               extra_search%log_bufs(index)%vector => lbuf

            endif
            !
            ! For collecting the data during the prism_put phase
            !
            send_info%nrecv                = send_info%nrecv + 1
            send_info%num2recv             = send_info%num2recv + rextra_msg(5,index)
            send_info%sender_global(index) = sender
            send_info%len_sent(index)      = rextra_msg(5,index)
            send_info%msg_id(index)        = rextra_msg(7,index)

         enddo ! loop over len pending Irecv requests

         !
         ! Add new cells to neighcells_3d array
         !
         len = sum(rextra_msg(3,:))

         Allocate(tgt_cell(len), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = len
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'tgt_cell, ... ', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         mm = 0
         nn = 0

         do n = 1, nrecv

            mm = mm + rextra_msg(5,n)

            do m = 1, rextra_msg(3,n)
               tgt_cell (nn+m) = tgt_ibuf(n)%vector(2*m)
            enddo

            nn = nn + rextra_msg(3,n)

         enddo

         len = sum(rextra_msg(5,:))

         Allocate(tgt_list(len), tgt_index(len), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 2*len
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'tgt_list, ... ', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         nn = 1

         do n = 1, nrecv
            mm = 1
            do m = 1, rextra_msg(3,n)
               do l = 1, tgt_ibuf(n)%vector(2*m)
                  tgt_list (nn) = tgt_ibuf(n)%vector(2*m-1)
                  tgt_index(nn) = nn
                  nn = nn + 1
                  mm = mm + 1
               enddo
            enddo
         enddo

         if ( nrecv > 1 ) &
              call psmile_quicksort_index ( tgt_list, len, tgt_index )

      endif ! nrecv > 0

      if ( nrecv > 0 .or. nsend > 0 ) then

         len = nbr_cells_tot+send_info%num2recv ! nbr_cells_tot+sum(rextra_msg(5,:))

         if ( len > 0 ) then

         Allocate(new_neighcells_3d(ndim_3d,len,1), stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = ndim_3d*len
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'new_neighcells_3d', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         l  = 1
         mm = 0
         nn = 0

         if ( nrecv > 0 ) then
             len = sum(rextra_msg(5,:))
         else
             len = 0
         endif

         do n = 1, ncpl

            do idx1 = 1, new_nbr_cells(n)
               new_neighcells_3d(1,mm+idx1,1) = neighcells_3d(1,nn+idx1,1)
               new_neighcells_3d(2,mm+idx1,1) = neighcells_3d(2,nn+idx1,1)
               new_neighcells_3d(3,mm+idx1,1) = neighcells_3d(3,nn+idx1,1)
            enddo

            nn = nn + nbr_cells(n)
            mm = mm + new_nbr_cells(n)

            if ( l <= len ) then

               do while ( n == tgt_list(l) )
                  !
                  ! as the extra cells are stored in a simple list at the
                  ! end, only one index is required, stored in place 2.
                  ! place 1 is used as identifier (extra_search%global_marker), place 3
                  ! is not used.
                  !
                  mm = mm + 1
                  new_nbr_cells(n) = new_nbr_cells(n) + 1
                  new_neighcells_3d(1,mm,1) = extra_search%global_marker
                  new_neighcells_3d(2,mm,1) = tgt_index(l)
                  new_neighcells_3d(3,mm,1) = PSMILe_Undef

                  l = l + 1

                  if ( l > len ) exit

               enddo ! while

            endif ! l <= len

         enddo ! n = 1, ncpl

         Deallocate(neighcells_3d, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = ndim_3d*nbr_cells_tot
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'neighcells_3d', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         neighcells_3d => new_neighcells_3d

         else

         Nullify(neighcells_3d)

         endif

      endif
      !
      ! Finally set nbr_cells to new values
      !
      nbr_cells = new_nbr_cells
!
!-------------------------------------------------------------------------
! 5.) Free memory
!-------------------------------------------------------------------------
!
      if ( max(nsend,nrecv) > 0 ) then

         ! In a loop over nsend the tgt_ibuf vector
         ! was already deallocated above.

         do n = 1, nrecv
            Deallocate (tgt_ibuf(n)%vector, stat=ierror)
            if (ierror /= 0) then
               ierrp (1) = PSMILe_Undef
               ierror = PRISM_Error_Dealloc
               call psmile_error ( ierror, 'tgt_ibuf', ierrp, 1, &
                    __FILE__, __LINE__ )
               return
            endif
         enddo

         Deallocate (tgt_ibuf, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) =  max(nsend,nrecv)
            ierror = PRISM_Error_Dealloc
            call psmile_error ( ierror, 'tgt_ibuf', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

      endif

      if ( nrecv > 0 ) then
         Deallocate(tgt_list, tgt_index, tgt_cell, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = 2*sum(rextra_msg(5,:)) + sum(rextra_msg(3,:))
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'tgt_list, ... ', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif

         Deallocate(rextra_msg, lrequest, stat=ierror)
         if (ierror /= 0) then
            ierrp (1) = (n_extra+1)*nrecv
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'rextra_msg, ... ', ierrp, 1, &
                 __FILE__, __LINE__ )
            return
         endif
      endif
!
!-------------------------------------------------------------------------
! 6.) All done
!-------------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

      call psmile_flushstd
#endif /* VERBOSE */

!
!  Formats:
!

#ifdef VERBOSE

9990 format (1x, a, ': psmile_global_search_cell_real: grid_id', i3, &
                    ' to ', i3, '(', i2, ')')
9980 format (1x, a, ': psmile_global_search_cell_real: eof ierror =', i3)

#endif /* VERBOSE */

      end subroutine PSMILe_global_search_cell_real