psmile_global_search_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_dble
!
! !INTERFACE:

      subroutine psmile_global_search_dble (comp_info,        &
                    control, len_cpl,                         &
                    var_id, grid_valid_shape,                 &
                    search, tgt_grid,                         &
                    neighbors_3d, nloc, num_neigh,            &
                    extra_search,                             &
                    interpolation_methods,                    &
                    interpolation_search, n_intmethods,       &
                    send_index,                               &
                    mask_available, use_mask, use_how,        &
                    grid_type, 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_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
!
      Type (dble_vector), Intent (In) :: tgt_grid (ndim_3d)

!     Coordinates to be searched

      Integer,            Intent (In) :: control (2, ndim_3d, search%search_data%npart)

!     Index range found for partition "ipart".

      Integer,            Intent (In) :: len_cpl (search%search_data%npart)

!     Number of points to be sent to the coupler for each partition.

      Integer,            Intent (In) :: var_id

!     Handle to the grid function

      Integer,            Intent (In) :: grid_valid_shape (2, ndim_3d)
      
!     Specifies the valid block shape for the "ndim_3d"-dimensional block

      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
!
      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) :: 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 ?

      Logical,            Intent (In) :: mask_available
!  !!!!! Fehlt da nicht noch ein "globales" mask_available !!!!
!
      Logical,            Intent (In) :: use_mask

      Integer,            Intent (In) :: use_how

!     Information about how to apply the mask

      Integer,            Intent (In) :: grid_type
!
!     Code for the type of grid generation
!     PRISM_Reglonlatvrt            = Regular grid in all directions
!     PRISM_Irrlonlat_regvrt        = Horizontal irregular, regular in the vertical
!     PRISM_Irrlonlatvrt            = Irregular grid in all directions
!     PRISM_Gaussreduced_regvrt     = Gaussian reduced grid with regular vertical
!
! !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_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS
!
!  NC_REG = Number of corners for regular directions
!  masked_out = Hash value for a point which is masked out
!
      Double Precision, Parameter  :: eps = 0.0d0
!
      Integer, Parameter           :: nc_reg = 2
      Integer, Parameter           :: masked_out = 0

!
! !LOCAL VARIABLES
!
!     ... Method pointer
!
      Integer                         :: comp_id, method_id
      Type (Method), Pointer          :: mp
!
      Type (Corner_Block), Pointer    :: corner_pointer
!
!     ... for the grid
!
      Integer                         :: grid_id, global_grid_id
      Type(Grid), Pointer             :: grid_info
!
      Type (integer_vector), Pointer  :: indices_req (:)
      Type (integer_vector), Pointer  :: required (:)
      Integer, Pointer                :: len_req (:)
!
!     ... Info's on actual component
!
      Integer                         :: rank, igrid

!
!     ... lookup table
!
      Integer, allocatable            :: igrid_to_dest_comp(:)
!
!     ... loop variables
!
      Integer                         :: ipart, j, n
!
!     ... for data on faces to other processes
!
      Logical                         :: mask_changed
      Integer                         :: nextra_prev, nreq, nx
!
      Integer, Allocatable            :: n_faces (:)
      Logical, Allocatable            :: send_mask (:)
      Double Precision, Allocatable   :: faces (:, :, :)
!
!     ... transformed faces
!
      Integer                         :: n_trans, n_int
!
      Integer                         :: tr_codes (max_num_trans_parts)
      Integer                         :: found (2*max_num_trans_parts)
      Real (PSMILe_float_kind)        :: face_range (2, ndim_3d)
      Real (PSMILe_float_kind)        :: dinter_trans (2, ndim_3d)
      Real (PSMILe_float_kind)        :: face_range_transformed (2, ndim_3d,           
                                                                 max_num_trans_parts), 
                                         dinter (2, ndim_3d, 2*max_num_trans_parts)
!
!     ... for points found
!
      Integer                         :: n_found, n_liste
!
!     ... for send to partners
!
      Type (Send_information), Pointer :: send_info
      Type (enddef_msg_extra)         :: extra_msg
      Integer                         :: extra_msg_buf (msg_extra_size)
      Integer                         :: answer (msg_extra_size)
      Logical                         :: virtual_cell_available
!
      Integer                         :: interp
      Integer                         :: i, ip, ipi, ipi_prev, ireq
      Integer                         :: len_ibuf, len_item, ndibuf, n_send
      Integer                         :: len_rbuf, len_rtem, ndrbuf, nd_len
      Integer, Pointer                :: len_nsend (:, :)
      Integer, Pointer                :: old_len (:, :)
      Integer, Allocatable            :: ibuf (:)
      Integer, Allocatable            :: srcloc_ind (:, :)
      Integer, Pointer                :: virtual_ind (:)
      Integer, Target                 :: dummy_virtual_ind (1)
      Double Precision, Pointer       :: buf (:)

!
!     ... for communication
!
      Integer                         :: dest, index, dest_comp, sender
      Integer                         :: nrecv, nprev_recv, 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_dble" performs a global search
! for the interpolation points not found in the local search.
!
! TODO: tgt_grid statt coords uebergeben;
!       kann man sich dann die unterschiedlichen Routine fuer
!       Get_faces_..., store_faces__ sparen !?!?
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 11.11.04    H. Ritzdorf  created
! 06.04.10    M. Hanke     adjusted in order to support multiple
!                          local grid id's for a single global grid id
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_global_search_dble.F90 3174 2011-05-04 09:41:43Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_global_search_dble.F90 3174 2011-05-04 09:41:43Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), var_id, search%msg_intersections%field_info%tgt_method_id, &
                               search%sender

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      ndrbuf = 0
      ndibuf = 0
      n_recv_req = 0
!
      method_id = Fields(var_id)%method_id
      mp => Methods(method_id)
!
      grid_id   = mp%grid_id
      grid_info => Grids(grid_id)
!
      global_grid_id = grid_info%global_grid_id
!
! Uebergeben oder ein-dimensional
!
      send_info => mp%send_infos_coupler(send_index)
      corner_pointer => grid_info%corner_pointer
!
      comp_id = comp_info%comp_id
!
      indices_req => extra_search%indices_req
      required    => extra_search%required
      len_req     => extra_search%len_req
      nreq        =  extra_search%n_req
!
      virtual_cell_available = Associated (send_info%virtual)
!
#ifdef PRISM_ASSERTION
      if (grid_type /= grid_info%grid_type) then
! nicht uebergeben !
         print *, 'grid_type', grid_type, grid_info%grid_type
         call psmile_assert ( __FILE__, __LINE__, &
                 'inconsistent grid types')
      endif
#endif
!
!  Rank of actual process in Component communicator
!
      rank    = Comps(comp_id)%rank
!
      if (count(interpolation_search(:)) > 1) then
         ierrp (1) = count(interpolation_search(:))

         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'Global search currently 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
!
! Can man die offset und partition information ausnutzen ?
!
!===> Allocate work vectors
!
!  n_faces = Number of block faces on which cell/control volume "i" is
!            located.
!  !!! n_faces wieder raus !!!
!  faces   = Range on block faces which interfere with faces of
!            cell/control volume "i".
!
      Allocate (faces (2, ndim_3d, nreq), STAT = ierror)

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

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'faces', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      Allocate (n_faces (nreq), STAT = ierror)

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

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

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

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'send_mask', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      mask_changed = .true.
!
!===> Get number of faces and intersection with block boundary
!     for each cell to be searched
!
!  len_item = Number of pieces of information
!             in integer buffer "ibuf"
!             ibuf(1:3) = srcloc of extra point
!                     (index of lower left corner)
!             ibuf(4)   = Index of the extra point to be searched
!                         This is a local index in the vectors
!                         "required" and "indices_req" of the actual
!                         partition "ipart".
!             ibuf(5)   = Code for interpolation bases required
!             ibuf(6)   = Virtual cell info (Only GaussReduced)
!             TODO: ibuf (4, *) raus
!  len_rtem = Number of pieces of information
!             per point in Double Precison buffer "rbuf"
!
! Ausserdem sollte man mehrfache srclocs rausfiltern
!
      nx = nreq ! oder weniger, wenn nicht alle
                ! cellen gesucht werden sollen
      len_item = ndim_3d + 2 ! default
!
      select case (grid_type)

         case (PRISM_Reglonlatvrt) 
            len_rtem = ndim_3d * (nc_reg + 1)
!
            call psmile_get_faces_3d_reg_dble (            &
                    search, extra_search,                  &
                    corner_pointer%corners_dble(1)%vector, &
                    corner_pointer%corners_dble(2)%vector, &
                    corner_pointer%corners_dble(3)%vector, &
                    corner_pointer%corner_shape,           &
                    corner_pointer%nbr_corners,            &
                    grid_valid_shape,                      &
                    neighbors_3d, nloc, num_neigh,         &
                    faces, n_faces, nreq, ierror)

         case (PRISM_Irrlonlat_regvrt) 
            len_rtem = ndim_3d + &
               (corner_pointer%nbr_corners/nc_reg)*ndim_2d + ndim_2d
!
            call psmile_get_faces_irreg2_dble (            &
                    search, extra_search,                  &
                    corner_pointer%corners_dble(1)%vector, &
                    corner_pointer%corners_dble(2)%vector, &
                    corner_pointer%corners_dble(3)%vector, &
                    corner_pointer%corner_shape,           &
                    corner_pointer%nbr_corners,            &
                    grid_valid_shape,                      &
                    neighbors_3d, nloc, num_neigh,         &
                    faces, n_faces, nreq, ierror)

         case (PRISM_Irrlonlatvrt) 
            len_rtem = ndim_3d + corner_pointer%nbr_corners*ndim_3d
!
            call psmile_get_faces_3d_dble (                &
                    search, extra_search,                  &
                    corner_pointer%corners_dble(1)%vector, &
                    corner_pointer%corners_dble(2)%vector, &
                    corner_pointer%corners_dble(3)%vector, &
                    corner_pointer%corner_shape,           &
                    corner_pointer%nbr_corners,            &
                    grid_valid_shape,                      &
                    neighbors_3d, nloc, num_neigh,         &
                    faces, n_faces, nreq, ierror)

         case (PRISM_Gaussreduced_regvrt) 
            if (virtual_cell_available) len_item = ndim_3d + 3
            len_rtem = ndim_3d * (nc_reg + 1)
!
            call psmile_get_faces_gauss2_dble (            &
                    search, extra_search,                  &
                    corner_pointer%corners_dble(1)%vector, &
                    corner_pointer%corners_dble(2)%vector, &
                    corner_pointer%corners_dble(3)%vector, &
                    corner_pointer%corner_shape,           &
                    corner_pointer%nbr_corners,            &
                    grid_valid_shape,                      &
                    neighbors_3d, nloc, num_neigh,         &
                    faces, nreq, ierror)

         case default
            ierrp (1) = grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'unsupported grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
            return

      end select
!
      if (ierror > 0) return
!
!===> Get range of faces and transform range into cyclic range
!
      face_range (1, 1) = minval (faces (1, 1, 1:nreq))
      face_range (2, 1) = maxval (faces (2, 1, 1:nreq))
      face_range (1, 2) = minval (faces (1, 2, 1:nreq))
      face_range (2, 2) = maxval (faces (2, 2, 1:nreq))
      face_range (1, 3) = minval (faces (1, 3, 1:nreq))
      face_range (2, 3) = maxval (faces (2, 3, 1:nreq))
!
      call psmile_transform_extent_cyclic (grid_info%grid_type, face_range,  &
                                           face_range_transformed, tr_codes, &
                                           n_trans, ierror)
      if (ierror > 0) return
!
!
!
      nd_len = 8
      Allocate (len_nsend (search%search_data%npart, nd_len), stat = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = search%search_data%npart * nd_len

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'len_nsend', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!     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
            send_mask (:) = .false.
            mask_changed = .false.
         endif
!
!        if current block is not part of the same global grid
         if (global_grid_id /= comp_info%all_extent_infos(igrid)%global_grid_id) cycle
!
!        if the current block is the one we are doing the global search for
!        (note ranks start with 0, ...)
         if (rank == dest_comp-1 .and. &
               comp_info%all_extent_infos(igrid)%local_grid_id == grid_id) cycle
!
!===> ... Does the "igrid"-th grid and the range of faces coincide ?
!
!  n_int = Number of significant intersections
!
!cdir vector
         do i = 1, n_trans

            dinter (1, :, i) = max (face_range_transformed(1,:,i), &
                                    comp_info%all_extent_infos(igrid)%extent(1,:))
            dinter (2, :, i) = min (face_range_transformed(2,:,i), &
                                    comp_info%all_extent_infos(igrid)%extent(2,:))

            ! in addition to the intersections we need to check whether the grid
            ! extent and the transformed face range has a common face at the
            ! common_grid_range boundaries (normally at -180/180)

            if (      face_range_transformed(1,1,i) == common_grid_range(1,1) .and. &
                .not. face_range_transformed(2,1,i) == common_grid_range(2,1)) then

               dinter (:, :, i+n_trans) = dinter (:, :, i)
               dinter (:, 1, i+n_trans) = common_grid_range(1,1)

            else if (.not. face_range_transformed(1,1,i) == common_grid_range(1,1) .and. &
                           face_range_transformed(2,1,i) == common_grid_range(2,1)) then

               dinter (:, :, i+n_trans) = dinter (:, :, i)
               dinter (:, 1, i+n_trans) = common_grid_range(2,1)

            else
               ! set dummy data
               dinter (1,:,i+n_trans) = 0
               dinter (2,:,i+n_trans) = -1
            endif
         end do ! i = 1, n_trans
!
         n_int = 0
!cdir vector
         do i = 1, 2 * 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 control volumes which
!         have intersection with the extent of the "n"-th grid
!???  What about cyclic boundaries ????
!???  use global offsets here !
!
         do i = 1, n_int
#ifdef DEBUGX
            print *, 'extent controlled', i, n_int, dinter(:,:,found(i))
#endif
!cdir vector
            do j = 1, nx
               send_mask (j) = send_mask (j) .or. &
                     (faces(2,1,j) >= dinter(1,1,found(i)) .and. &
                      faces(1,1,j) <= dinter(2,1,found(i)) .and. &
                      faces(2,2,j) >= dinter(1,2,found(i)) .and. &
                      faces(1,2,j) <= dinter(2,2,found(i)) .and. &
                      faces(2,3,j) >= dinter(1,3,found(i)) .and. &
                      faces(1,3,j) <= dinter(2,3,found(i)))
#ifdef DEBUGX
print *, j, faces(:,:,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
         if (n_recv_req > nd_len) then
            nd_len = nd_len + 16

!           Reallocate len_nsend

            old_len => len_nsend
            Allocate (len_nsend (search%search_data%npart, nd_len), stat = ierror)

            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = search%search_data%npart * nd_len

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'len_nsend', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
!
            len_nsend (:, 1:n_recv_req-1) = old_len (:, 1:n_recv_req-1)
!
            Deallocate (old_len)
         endif
!
#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
!
!  Double precision vector buf (1:ndrbuf):
!  buf (:) = Coordinates of the point to be searched
!            ndim_3d data items
!  buf (:) = Coordinates of corners of control volume/cell
!            "len_rtem" data items
!
!  Integer buffer ibuf (1:ndibuf):
!
!  ibuf (:) = srcloc of extra point
!             ndim_3d data items
!  ibuf (:) = Index of the extra point to be searched
!             1       data items
!
         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
!
!===> Get Indices (in srclocs) of points to be sent
!     TODO: Direkt in ibuf speichern !
!
         Allocate (srcloc_ind(ndim_3d, n_send), STAT = ierror)

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

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'srcloc_ind', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
         select case (send_info%nvec)
         case (1)
            call psmile_get_face_ind_3d  (search, extra_search,         &
                           send_info, len_cpl,                          &
                           send_mask, nreq, srcloc_ind, n_send,         &
                           ierror)

         case (2)
            call psmile_get_face_ind_21d (search, extra_search,         &
                           send_info, len_cpl,                          &
                           send_mask, nreq, srcloc_ind, n_send,         &
                           ierror)
!
         case (ndim_3d)
            call psmile_get_face_ind_reg (search, extra_search,         &
                           send_info, len_cpl,                          &
                           send_mask, nreq, srcloc_ind, n_send,         &
                           ierror)
         end select
!
         if (ierror > 0) return
!
         if (grid_type == PRISM_Gaussreduced_regvrt) then
#ifdef DEBUG
! das wurde fuer die alten Gaussreduced Gird (vor virtual_cell info)
! gebraucht; derzeit noch drin wg. conservative
  print *, '### psmile_global_search_dble.F90 raus damit !!',  send_info%nvec
#endif
            do i = 1, n_send
               srcloc_ind (1, i) = abs (srcloc_ind(1,i))
            end do
         endif
!
         if (virtual_cell_available) then
            Allocate (virtual_ind(n_send), STAT = ierror)

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

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'virtual_ind', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
!
            call psmile_get_faces_virtual_ind (search, extra_search, &
                      send_info, len_cpl,                            &
                      send_mask, nreq, virtual_ind, n_send,          &
                      ierror)
            if (ierror > 0) return

         else

            virtual_ind => dummy_virtual_ind

         endif
!
!===> Store data in send buffers
!     ??? compact machen !?!?
!         Mehrfache Zellen Info's in ibuf und corner volumen
!         eliminieren ?
!     TODO: ipi und ip wieder vereinigen
!
         nextra_prev = 0
         ipi = 0
         ip  = 0
!
         select case (grid_type)

         case (PRISM_Reglonlatvrt) 
            do ipart = 1, search%search_data%npart
!
            ipi_prev = ipi
!
            if (len_req (ipart) > 0) then
!
               call psmile_store_faces_3d_reg_dble (           &
                        indices_req(ipart)%vector,             &
                        required   (ipart)%vector,             &
                        len_req (ipart),                       &
                        tgt_grid (1)%vector,                   &
                        tgt_grid (2)%vector,                   &
                        tgt_grid (3)%vector, nloc,             &
                        corner_pointer%corners_dble(1)%vector, &
                        corner_pointer%corners_dble(2)%vector, &
                        corner_pointer%corners_dble(3)%vector, &
                        corner_pointer%corner_shape,           &
                        corner_pointer%nbr_corners,            &
                        grid_valid_shape,                      &
                        send_mask(nextra_prev+1:),             &
                        srcloc_ind,                            &
                        ibuf, len_item, n_send, ipi,           &
                        buf,  len_rtem, n_send, ip, ierror)

                nextra_prev = nextra_prev + len_req (ipart)
             endif
!
             len_nsend (ipart, n_recv_req) = ipi - ipi_prev
             end do
!
#ifdef PRISM_ASSERTION
         if (ipi /= n_send) then
            print *, 'ipi, n_send', ipi, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer ibuf')
         endif
!
         if (ip /= n_send) then
            print *, 'ip, n_send', ip, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer buf')
         endif
#endif
!
         case (PRISM_Irrlonlat_regvrt) 
!
            do ipart = 1, search%search_data%npart
!
            ipi_prev = ipi
!
            if (len_req (ipart) > 0) then
!
               call psmile_store_faces_irreg2_dble (           &
                        indices_req(ipart)%vector,             &
                        required   (ipart)%vector,             &
                        len_req (ipart),                       &
                        tgt_grid (1)%vector,                   &
                        tgt_grid (2)%vector,                   &
                        tgt_grid (3)%vector, nloc,             &
                        corner_pointer%corners_dble(1)%vector, &
                        corner_pointer%corners_dble(2)%vector, &
                        corner_pointer%corners_dble(3)%vector, &
                        corner_pointer%corner_shape,           &
                        corner_pointer%nbr_corners,            &
                        grid_valid_shape,                      &
                        send_mask(nextra_prev+1:),             &
                        srcloc_ind,                            &
                        ibuf, len_item, n_send, ipi,           &
                        buf,  len_rtem, n_send, ip, ierror)

                nextra_prev = nextra_prev + len_req (ipart)
             endif
!
             len_nsend (ipart, n_recv_req) = ipi - ipi_prev
             end do
!
#ifdef PRISM_ASSERTION
         if (ipi /= n_send) then
            print *, 'ipi, n_send', ipi, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer ibuf')
         endif
!
         if (ip /= n_send) then
            print *, 'ip, n_send', ip, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer buf')
         endif
#endif
!
         case (PRISM_Gaussreduced_regvrt) 
!
            do ipart = 1, search%search_data%npart
!
            ipi_prev = ipi
!
            if (len_req (ipart) > 0) then
!
               call psmile_store_faces_gauss2_dble (           &
                        indices_req(ipart)%vector,             &
                        required   (ipart)%vector,             &
                        len_req (ipart),                       &
                        tgt_grid (1)%vector,                   &
                        tgt_grid (2)%vector,                   &
                        tgt_grid (3)%vector, nloc,             &
                        corner_pointer%corners_dble(1)%vector, &
                        corner_pointer%corners_dble(2)%vector, &
                        corner_pointer%corners_dble(3)%vector, &
                        corner_pointer%corner_shape,           &
                        corner_pointer%nbr_corners,            &
                        grid_id, grid_valid_shape,             &
                        send_mask(nextra_prev+1:),             &
                        srcloc_ind, virtual_ind,               &
                        virtual_cell_available,                &
                        ibuf, len_item, n_send, ipi,           &
                        buf,  len_rtem, n_send, ip, ierror)

                nextra_prev = nextra_prev + len_req (ipart)
             endif
!
             len_nsend (ipart, n_recv_req) = ipi - ipi_prev
             end do
!
#ifdef PRISM_ASSERTION
         if (ipi /= n_send) then
            print *, 'ipi, n_send', ipi, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer ibuf')
         endif
!
         if (ip /= n_send) then
            print *, 'ip, n_send', ip, n_send
            call psmile_assert ( __FILE__, __LINE__, &
                    'inconsistent length for pack buffer buf')
         endif
#endif

         case (PRISM_Irrlonlatvrt) 
            do ipart = 1, search%search_data%npart
!
            ipi_prev = ipi
!
            if (len_req (ipart) > 0) then
!
               call psmile_store_faces_3d_dble (               &
                        indices_req(ipart)%vector,             &
                        required   (ipart)%vector,             &
                        len_req (ipart),                       &
                        tgt_grid (1)%vector,                   &
                        tgt_grid (2)%vector,                   &
                        tgt_grid (3)%vector, nloc,             &
                        corner_pointer%corners_dble(1)%vector, &
                        corner_pointer%corners_dble(2)%vector, &
                        corner_pointer%corners_dble(3)%vector, &
                        corner_pointer%corner_shape,           &
                        corner_pointer%nbr_corners,            &
                        grid_valid_shape,                      &
                        send_mask(nextra_prev+1:),             &
                        srcloc_ind,                            &
                        ibuf, len_item, n_send, ipi,           &
                        buf,  len_rtem, n_send, ip, ierror)

                nextra_prev = nextra_prev + len_req (ipart)
             endif
!
             len_nsend (ipart, n_recv_req) = ipi - ipi_prev
             end do

         case default
            ierrp (1) = grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'unsupported grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
            return

         end select
!
         Deallocate (srcloc_ind)
         if (virtual_cell_available) Deallocate (virtual_ind)
!
#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 = extra_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.
!
         call psmile_init_enddef_msg_extra(extra_msg)

         extra_msg%reqest_type         = interpolation_methods (interp)
         extra_msg%datatype            = PRISM_DOUBLE_PRECISION
         extra_msg%len_int_data        = len_ibuf
         extra_msg%len_coord_data      = len_rbuf
         extra_msg%global_comp_id      = comp_info%global_comp_id
         extra_msg%transi_out_id       = search%msg_intersections%field_info%transient_out_id
!
         extra_msg%num_volumes         = n_send
         extra_msg%num_int_per_vol     = len_item
         extra_msg%num_items_per_coord = len_rtem
!
         extra_msg%idx_req             = n_recv_req
         extra_msg%num_neigh           = num_neigh

         extra_msg%local_grid_id = comp_info%all_extent_infos(igrid)%local_grid_id
!
#ifdef DEBUGX
         print *, 'grid shape', grid_valid_shape
         do i = 1, n_send
            print *, ibuf ((i-1)*len_item+1:i*len_item)
         end do
#endif
!
         if (Associated (grid_info%partition)) then
!           First dimension: Blocks
!           Last  dimension: Global indices of partition
!
            extra_msg%partition_avail = .true.
            extra_msg%partition(:)    = grid_info%partition (1, 1:3)
!
!===> ... Transfer local coordinates into global coordinates
!
            if (grid_info%grid_type == PRISM_Gaussreduced_regvrt) then

               call psmile_trans_loc2glob_gauss2 (grid_id, &
                         ibuf, len_item, n_send, ierror)

            else
!
               call psmile_trans_loc2glob_3d (grid_id,     &
                         ibuf, len_item, n_send, ierror)

            endif
!
            if (ierror > 0) return

         else
!
!           ... No partition info
!
            extra_msg%partition_avail = .false.
         endif
!
#ifdef DEBUG
         print *, ' Sending ', exttag, ' to destination ', dest
#endif
         call psmile_pack_msg_extra (extra_msg, extra_msg_buf)

         call psmile_bsend (extra_msg_buf, msg_extra_size, 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
!
#ifdef DEBUG
         print *, ' Sending ', exttag, ' to destination ', dest
#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
!
#ifdef DEBUG
         print *, ' Sending ', exttag, ' to destination ', dest
#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
            ! MoHa: Be carefull even though this message uses rexttag and the size
            ! msg_extra_size it does not contain the same information as
            ! enddef_msg_extra! This messages is sent by psmile_return_extra_off.
            ! There is no documentation on the exact content of this message!
            ! (From my point of view, this is "no good" programming style.)
            call MPI_Irecv (answer, msg_extra_size, 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, n_faces, faces, igrid_to_dest_comp)
!
      if (ndrbuf > 0) Deallocate (buf)
      if (ndibuf > 0) Deallocate (ibuf)
!
!-------------------------------------------------------------------------
!     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
! nprev_recv = Number of data items (per coordinate) received
!              in prevoius receives.
!
      nprev_recv = 0
      nrecv = 0
!
      if (n_recv_req > 0) then
         Allocate (extra_search%dble_bufs(n_recv_req),  &
                   send_info%sender_global(n_recv_req), &
                   send_info%len_sent(n_recv_req),      &
                   send_info%msg_id(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, 'extra_search%dble_bufs', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif
      endif
!
      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), ': psmile_global_search_dble: action%nreq, recv_req', &
                     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) /= interpolation_methods (interp)) then
                print *, 'answer(1)', answer (1), interpolation_methods (interp)
               call psmile_assert (__FILE__, __LINE__, &
                    "incorrect message received ? interpolations doesn't fit")
            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
!
            extra_search%dble_bufs (nrecv)%vector  => buf
            send_info%sender_global (nrecv) = sender
            send_info%len_sent      (nrecv) = n_liste
            send_info%msg_id        (nrecv) = answer (16)
!
            if (nprev_recv > 0) then
!              Increment index_found (for collected list of buf's)
!              except for point masked out
!cdir vector
               do i = 2*n_send+1, 2*n_send+n_found
                  if (ibuf(i) /= masked_out) then
                     ibuf (i) = ibuf (i) + nprev_recv
                  endif
               end do
            endif
!
            nprev_recv = nprev_recv + n_liste
!
!===> ... Add information on points found by global search 
!
            call psmile_add_points_found (grid_id, search, extra_search, &
                        ibuf (1:n_send), ibuf (n_send+1:2*n_send),       &
                        n_send, len_nsend(1:search%search_data%npart,ireq),          &
                        ibuf (2*n_send+1:2*n_send+n_found), n_found,     &
                        neighbors_3d, nloc, num_neigh,                   &
                        grid_valid_shape, use_how, ierror)
                   if (ierror > 0) return
!
            Deallocate (ibuf)
         endif
!
!===> ... Setup new request
!
         if (n < n_recv_req) then
            ! MoHa: Be carefull even though this message uses rexttag and the size
            ! msg_extra_size it does not contain the same information as
            ! enddef_msg_extra! This messages is sent by psmile_return_extra_off.
            ! There is no documentation on the exact content of this message!
            ! (From my point of view, this is "no good" programming style.)
            call MPI_Irecv (answer, msg_extra_size, 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
!
!    ... 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)
!
      Deallocate (len_nsend)
!
!===> Save number of data to be received in "extra_search"
!
      send_info%nrecv    = nrecv
      send_info%num2recv = nprev_recv
!
#ifdef PRISM_ASSERTION
      if (nrecv > 0) then
         if (send_info%num2recv /= SUM (send_info%len_sent(1:nrecv))) then
            print *, 'num2recv', send_info%num2recv, &
                     SUM (send_info%len_sent(1:nrecv)), nrecv
            call psmile_assert (__FILE__, __LINE__, &
                 "inconsistent num2recv")
         endif
      endif
#endif
!
!-------------------------------------------------------------------------
!===> 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_dble: var_id', i3, &
                    ' to ', i3, '(', i2, ')')
9980 format (1x, a, ': psmile_global_search_dble: eof ierror =', i3)
9970 format (1x, a, ': psmile_global_search_dble: send from', i3, &
                    ' to', i3, '[', i3, '], n_send =', i6)

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

#ifdef DEBUG
#endif

      end subroutine PSMILe_global_search_dble