psmile_info_trs_loc_irreg2_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_info_trs_loc_irreg2_dble
!
! !INTERFACE:

      subroutine psmile_info_trs_loc_irreg2_dble (comp_info,  &
                           coords, shape, control, len_cpl,   &
                           var_id, grid_valid_shape,          &
                           search, method_id,                 &
                           send_index, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_info_trs_loc_irreg2_dble
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      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) :: coords (ndim_3d, search%npart)

!     Coordinates to be searched

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

!     Dimension of coordinate arrays "coords (1:ndim_3d, ipart)", 
!     which contain the coordinates to be searched for partition "ipart".
!     The dimension depends on search%grid_type !

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

!     Index range found for partition "ipart".

      Integer, Intent (InOut)         :: len_cpl (search%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)            :: method_id

!     Method id

      Integer, Intent (InOut)         :: send_index

!     Send index of coupler send for the method and coordinates

! !OUTPUT PARAMETERS:

      integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILe_info_trs_loc_irreg2_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
! indl = Index of LonLat values in "srcloc"
! indz = Index of Vert   values in "srcloc"
!
      Integer, Parameter              :: indl = 1
      Integer, Parameter              :: indz = 2
!
! !LOCAL VARIABLES
!
!     ... field pointer
!
      Type (Gridfunction), Pointer    :: field
!
!     ... Method pointer
!
      Type (Method), Pointer          :: mp
!
!     ... Grid pointer
!
      Type (Grid), Pointer            :: grid_info
!
      Integer                         :: grid_id
      Integer                         :: ncpl, new_ncpl
      Integer                         :: jj, totsiz, offset
!
!     ... loop variables
!
      Integer                         :: i, j, n, n1, n2, n3, nn
      Integer                         :: jpart
!
!     ... for search of interpolation bases
!
      Type (Corner_Block), Pointer    :: cp
      Type (Send_information), Pointer :: send_info
      Type (Extra_search_info)        :: extra_search
!
      Integer                         :: nb_neighbors, nb_extra
      Integer, Pointer                :: dstijk (:, :)
      Integer, Allocatable            :: new_dstijk (:, :)
!
      Integer, Allocatable            :: neighbors_3d (:, :, :)
      Integer                         :: interpolation(ndim_3d)
      Integer                         :: interpolation_mode, search_mode
      Integer                         :: interpolation_type
      Integer                         :: interpolation_methods (indz)
      Integer                         :: n_intmethods
      Logical                         :: interpolation_search  (indz)
      Logical                         :: global_search
      Logical                         :: global_cell_search
!
!     ... for masks
!         dummy_mask_array is a dummy array for to be
!         transferred to interpolation routines
!         Note: The target attributes (see !rr) were removed
!               since there problems with a compiler and array bound checking
!
      Integer                         :: id_src_mask, id_tgt_mask
      Integer                         :: mask_id
      Integer                         :: use_how (3)
      Logical                         :: src_mask_available, use_mask
!
      Logical, Target                 :: dummy_mask_array (1)
!rr   Integer, Target                 :: dummy_mask_shape (2, ndim_3d)
      Integer                         :: dummy_mask_shape (2, ndim_3d)
      Integer, Target                 :: dummy_mask (1) ! dummy vector
!
!rr   Integer, Pointer                :: mask_shape (:, :)
      Integer                         :: mask_shape (2, ndim_3d)
      Logical, Pointer                :: mask_array (:)
      Integer, Pointer                :: src_mask (:)
!
!     ... for partitions
!
      Integer                         :: ibeg, iend, ipart, nprev
      Integer                         :: new_len_cpl(search%npart)
!
!     ... for transfer to the coupler (source grid)
!
      Integer                         :: cpl_id
      Integer                         :: epio_id, trs_rank
      Integer                         :: local_src_size, shift, src_size
      Integer                         :: src_cell_size
      Integer                         :: id_trans_out
      Integer                         :: ip, len, leni, lenij
!
      Integer, Pointer                :: list_entries (:, :)
      Type (dble_vector)              :: src_grid (ndim_3d)
      Type (dble_vector)              :: src_cell (ndim_3d)
      Integer, Allocatable            :: neighbors  (:, :)
      Integer, Pointer                :: neighcells (:, :)
      Integer, Allocatable            :: source_cell_index (:)
      Integer, Pointer                :: nbr_cells (:)
      Integer, Pointer                :: new_nbr_cells (:)
      Integer                         :: nbr_cells_tot, nb_corners
      Integer                         :: nbr_corners(ndim_3d)
      Integer                         :: corner_shape    (2, ndim_3d)
      Integer                         :: corner_shape_3d (2, ndim_3d, ndim_3d)
!
!     ... for transfer to the coupler (target grid)
!
      Integer                         :: id_trans_in, dest_comp_id
      Logical                         :: use_target_grid
      Logical                         :: use_target_cell
      Logical                         :: use_source_cell
      Integer                         :: tgt_corners(3)
!
      Type (dble_vector)              :: tgt_grid (ndim_3d)
      Type (dble_vector)              :: tgt_cell (ndim_3d)
      Type (dble_vector)              :: new_tgt_cell (ndim_3d)
      Integer, Pointer                :: tgt_mask(:)
!
!     ... for transformed values
!
      logical                         :: transformed, cubic
      Type (dble_vector)              :: sinvec, cosvec
!
!     ... for error handling
!
      Integer, parameter              :: nerrp = 3
      Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_info_trs_loc_irreg2_dble" informs the coupler/
! transformer on the data to be interpolated
! for the method (grid) and the subgrid coords of the sending process
! "search%sender".
!
! Subroutine "PSMILe_info_trs_loc_irreg2_dble" is designed for grids of
! of type "PRISM_Irrlonlat_regvrt".
!
! !TO DO!
!
! 1.) neighbors and neighcells can be united
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_info_trs_loc_irreg2_dble.F90 2936 2011-02-03 09:36:47Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_info_trs_loc_irreg2_dble.F90 2936 2011-02-03 09:36:47Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization of dummy mask shape
!
      data dummy_mask_shape /0, 1, 0, 1, 0, 1/
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), method_id, search%msg_intersections%first_tgt_method_id, &
                               search%sender

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror  = 0
!
      field => Fields(var_id) 
      mp    => Methods(method_id)
      send_info => mp%send_infos_coupler(send_index)
      grid_id = mp%grid_id
      grid_info => Grids(grid_id)
      cp    => Grids(grid_id)%corner_pointer
!
      transformed = .false.
      cubic = .false.
!
      global_search = .false.
      global_cell_search = .false.
      interpolation_search (:) = .false.
!
      use_how(1) = PSMILe_undef
      use_how(2) = PSMILe_undef
      use_how(3) = grid_id
!
#ifdef PRISM_ASSERTION
!
!===> Control Assertions
!
      if ( .not. associated(field%Taskout) ) then
         print *, 'var_id', var_id
         call psmile_assert (__FILE__, &
                 __LINE__, "Taskout is not available for Field")
      end if

      if (var_id < 1 .or. var_id > Number_of_Fields_allocated) then
         print *, 'var_id', var_id
         call psmile_assert (__FILE__, &
                 __LINE__, "Incorrect var_id")
      endif

      if ( mp%status == PSMILe_Status_free .or. &
           mp%status == PSMILe_Status_undefined ) then
         call psmile_assert (__FILE__, &
                 __LINE__, "Incorrect method")
      endif

      if (send_index < 1) then
         print *, send_index
         call psmile_assert (__FILE__, &
                 __LINE__, "send_index is out of range")
      endif

! This assertion is not yet save for conservative remapping.
! Here it can happen that during psmile_store_dest_locs_3d
! all target cells which shall be added are masked out thus
! making nloc to zero. On the other hand hand, even for nloc
! zero we have to go through for conservative remapping to
! avoid dead locks. 
!
!      if (send_info%nloc < 1) then
!         print *, "ncpl", send_info%nloc
!         call PSMILe_Assert (__FILE__, &
!                 __LINE__, "ncpl > 0 expected")
!      endif

      if (send_info%nvec == indz) then
         if (search%grid_type /= PRISM_Reglonlatvrt     .and. &
             search%grid_type /= PRISM_Irrlonlat_Regvrt .and. &
             search%grid_type /= PRISM_Gaussreduced_Regvrt ) then
            print *, "nvec", send_info%nvec
            call psmile_assert (__FILE__, &
                    __LINE__, "nvec == indz not expected for this grid type")
         endif
      else if (send_info%nvec /= 1) then
         print *, "nvec", send_info%nvec
         call psmile_assert (__FILE__, &
                  __LINE__, "nvec must be in range [1:2]")
      endif

      if (send_info%nloc /= Sum(len_cpl)) then
         print *, "ncpl", send_info%nloc, Sum (len_cpl)
         call psmile_assert (__FILE__, &
                 __LINE__, "ncpl == SUM(len_cpl) expected")
      endif
#endif
!
!===> Get send information
!
      ncpl   =  send_info%nloc
      dstijk => send_info%dstijk
!
!===> Get information on donor (source) grid
!     use_mask = Should "mask_array" be used to compute nearest neighbours ?
!                .true.  : Use values in "mask_array" for search.
!                .false. : Corresponds to used_masked in XML file.
!                          Send mask to transformer.
!
      mask_id = field%mask_id
      src_mask_available = mask_id /= PRISM_UNDEFINED
      use_mask = src_mask_available
!
      if (src_mask_available) then
         mask_array => Masks(mask_id)%mask_array
!rr      mask_shape => Masks(mask_id)%mask_shape
         mask_shape =  Masks(mask_id)%mask_shape
      else
         mask_array => dummy_mask_array
!rr      mask_shape => dummy_mask_shape
         mask_shape =  dummy_mask_shape
      endif
!
!----------------------------------------------------------------------------
!     Get interpolation type requested
!----------------------------------------------------------------------------
!
      id_trans_in  = search%msg_intersections%transient_in_id
      id_trans_out = search%msg_intersections%transient_out_id
      dest_comp_id = all_comp_infos (search%msg_intersections%all_comp_infos_comp_idx)%global_comp_id
!
         do i = 1, size (field%Taskout)
         if (field%Taskout(i)%origin_type      == PSMILe_comp  .and. &
             field%Taskout(i)%remote_comp_id   == dest_comp_id .and. &
             field%Taskout(i)%remote_transi_id == id_trans_in  .and. &
             field%Taskout(i)%global_transi_id == id_trans_out) exit
         end do
!
      if (i > size (field%Taskout)) then
#ifdef DEBUG
         print *, 'psmile_info_trs_loc_irreg2_dble: method_id ', method_id
         print *, 'search for Taskout id ', id_trans_out, ' and Taskin Id ', &
                  id_trans_in, ' for dest component', dest_comp_id
!
            do i = 1, size (field%Taskout)
            print 9960, i, &
                        field%Taskout(i)%origin_type, &
                        field%Taskout(i)%remote_comp_id, &
                        field%Taskout(i)%remote_transi_id, &
                        field%Taskout(i)%global_transi_id
            end do
#endif /* DEBUG */
!
         ierrp (1) = id_trans_in
         ierrp (2) = id_trans_out
         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'cannot find Taskout index for transient ids', &
                             ierrp, 2, __FILE__, __LINE__ )
      endif
!
      interpolation_type    = field%Taskout(i)%interp%interp_type
      interpolation_methods = field%Taskout(i)%interp%interp_meth (1:2)
!
!===> Get interpolation type and interpolation method
!     interp_meth(1)   = 3d       Interpolation
!     interp_meth(1:2) = 2d,1d    Interpolation
!     interp_meth(1:3) = 1d,1d,1d Interpolation
!
#ifdef DEBUG
      print *, "interpolation type  ", interpolation_type
      print *, "interpolation method", interpolation_methods
#endif
!
      select case (interpolation_type)
         case (PSMILe_3D)
            n_intmethods = 1

         case (PSMILe_2D1D)
            n_intmethods = 2

!     Error: unsupported interpolation type

         case default
            ierrp (1) = interpolation_type
            ierror = PRISM_Error_Interp_type

            call psmile_error ( ierror, &
                'unsupported interpolation for gridtype PRISM_Irrlonlat_regvrt', &
                 ierrp, 1, __FILE__, __LINE__ )
            return

      end select
!
         do j = n_intmethods, 1, -1
         if (interpolation_methods(j) /= PSMILe_none) exit
         end do
!
      n_intmethods = j
!
      if (n_intmethods == 0) then
         ierror = PRISM_Error_Interp_type
         call psmile_error ( ierror, 'All interpolations are PSMILE_none', &
                             ierrp, 0, __FILE__, __LINE__ )
         return
      endif
!
!===> ... Get interpolation methods and number of neighbours
!
      nb_neighbors = 1
      nb_extra     = 1 ! number of nearest neighbours for extra search
!
      do j = 1, n_intmethods
         select case (interpolation_methods(j))
!
!           Conservative remapping
!
            case (PSMILe_conserv2D)
               nb_neighbors = nb_neighbors * 4
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               use_how(2)   = PSMILe_conserv2D
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global

            case (PSMILe_conserv3D)
               nb_neighbors = nb_neighbors * 8
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               use_how(2)   = PSMILe_conserv3D
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global
!
!           Nearest neighbour
!
            case (PSMILe_nnghbr3D, PSMILe_nnghbr2D)
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               nb_neighbors = nb_neighbors * field%Taskout(i)%interp%arg2(j)
               nb_extra     = nb_extra     * field%Taskout(i)%interp%arg2(j)
               use_mask     = use_mask .and. &
                              (field%Taskout(i)%interp%arg5(j) /= 1)
!
!           Linear search
!
            case (PSMILe_trilinear)
               nb_neighbors = nb_neighbors * 8
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global

            case (PSMILe_bilinear)
               nb_neighbors = nb_neighbors * 4
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global

            case (PSMILe_linear)
! Caution: This will overwrite the value for bilinear in case
!          bilinear & linear are active for a field
               nb_neighbors = nb_neighbors * 2
               use_how(1)   = field%Taskout(i)%interp%arg4(j)
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global
!
!           Cubic search
!
            case (PSMILe_bicubic)
               nb_neighbors = nb_neighbors * 4 * 4
               use_how(1)  = field%Taskout(i)%interp%arg4(j)
               cubic = .true.
               interpolation_search (j) = &
                  field%Taskout(i)%interp%arg3 (j) == PSMILE_global
!
            case (PSMILe_none)

            case (PSMILe_undef)
               ierrp (1) = j
               ierror = PRISM_Error_Interp_type

               call psmile_error ( ierror, &
                  'undefined interpolation method (PSMILE_UNDEF) found', &
                  ierrp, 1, __FILE__, __LINE__ )
               return

            case default
               ierrp (1) = interpolation_methods(j)
               ierror = PRISM_Error_Internal

               call psmile_error ( ierror, &
                  'unsupported interpolation for gridtype PRISM_Irrlonlat_regvrt', &
                  ierrp, 1, __FILE__, __LINE__ )
               return

         end select
      end do
!
      if (nb_neighbors < 1) then
         ierror = PRISM_Error_Interp_type
         ierrp (1) = nb_neighbors
!
         call psmile_error ( ierror, &
                 'Number of interpolation bases is less than 1', &
                 ierrp, 1, __FILE__, __LINE__ )
         return
      endif
!
      if (comp_info%size > 1) then
         if ( interpolation_methods(1) == PSMILe_conserv2D ) then
            global_cell_search = ANY (interpolation_search)
         else
            global_search = ANY (interpolation_search)
         endif
!
#ifdef DEBUG
         print *, "interpolation search", interpolation_search (:)
         print *, "global        search", global_search
         print *, "global cell   search", global_cell_search
#endif
      else
         interpolation_search = .false.
      endif
!
#ifdef PRISM_ASSERTION
         do j = 1, n_intmethods
         select case (interpolation_methods(j))
!
!           Conservative remapping
!           not yet supported
            case (PSMILe_conserv2D)
               if (n_intmethods > 2 .or. &
                  (n_intmethods == 2 .and. &
                   interpolation_methods(2) /= PSMILe_linear)) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_conserv2D with that interpolation")
               endif

            case (PSMILe_conserv3D)
               if (n_intmethods > 1) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_trilinear with other interpolations")
               endif
!
!           Nearest neighbour
!
            case (PSMILe_nnghbr3D)
               if (n_intmethods > 1) then
                  call psmile_assert (__FILE__, __LINE__, &
              "Don't know to combine PSMILe_nnghbr3D with other interpolations")
               endif

            case (PSMILe_nnghbr2D)
               if (n_intmethods > 1) then
                  call psmile_assert (__FILE__, __LINE__, &
              "Don't know to combine PSMILe_nnghbr2D with other interpolations")
               endif
!
!           Tri-linear interpolation
!
            case (PSMILe_trilinear)
               if (n_intmethods > 1) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_trilinear with other interpolations")
               endif
!
!           Bi-linear interpolation
            case (PSMILe_bilinear)
               if (n_intmethods > 2 .or. &
                  (n_intmethods == 2 .and. &
                   interpolation_methods(2) /= PSMILe_linear)) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_bilinear with that interpolation")
               endif
!
!           Linear interpolation
            case (PSMILe_linear)
               if (n_intmethods  < 2 .or.  &
                  (n_intmethods == 2 .and. &
                  (interpolation_methods(1) /= PSMILe_bilinear .and. &
                   interpolation_methods(1) /= PSMILe_bicubic))) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_linear with that interpolation")
               endif
!
!           Bi-cubic interpolation
!           Supported: Bicubic/Linear Interpolation
!
            case (PSMILe_bicubic)
               if (n_intmethods  > 2 .or.  &
                  (n_intmethods == 2 .and. &
                   interpolation_methods(2) /= PSMILe_linear)) then
                  call psmile_assert (__FILE__, __LINE__, &
               "Don't know to combine PSMILe_bicubic with that interpolation")
               endif

         end select
         end do
#endif

!----------------------------------------------------------------------------
!     Get coordinates of the target grid.
!     The coordinates are stored as single array for each direction.
!
! TODO: Here we have to control the Grid Types which can be transferred.
!
!----------------------------------------------------------------------------
!
      use_target_cell = ( interpolation_methods(1) == PSMILe_Conserv2D .or. &
                          interpolation_methods(1) == PSMILe_Conserv3D )

      use_source_cell = use_target_cell

      if (use_target_cell) then
!
!===> Conservative unterpolation using Cell data -----------------------
!
         use_target_grid = .false.

         if ( search%grid_type == PRISM_Irrlonlat_regvrt ) then
            tgt_corners(1:2) = 4
            tgt_corners(3)   = 4

            do i = 1, ndim_3d
               Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = ncpl

                  ierror = PRISM_Error_Alloc
                  call psmile_error ( ierror, 'tgt_cell(i)%vector', &
                       ierrp, 2, __FILE__, __LINE__ )
                  return
               endif
            end do
         else if ( search%grid_type == PRISM_Reglonlatvrt ) then

            tgt_corners(:) = 4

            do i = 1, ndim_3d
               Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = ncpl

                  ierror = PRISM_Error_Alloc
                  call psmile_error ( ierror, 'tgt_cell(i)%vector', &
                       ierrp, 2, __FILE__, __LINE__ )
                  return
               endif
            end do

         else if ( search%grid_type == PRISM_Gaussreduced_regvrt ) then

            tgt_corners = nb_neighbors

            nprev = 0
            ncpl  = 0

            do ipart = 1, search%npart

               totsiz =  size(coords(1,ipart)%vector)
               offset =  totsiz/2

               if (len_cpl (ipart) > 0) then
                  ibeg = nprev + 1
                  iend = nprev + len_cpl (ipart)

                  do j = ibeg, iend
                     if ( dstijk(1,j) /= PSMILe_Undef .and. &
                          dstijk(1,j)-shape(1,1,ipart)+1 <= offset ) ncpl = ncpl + 1

!rr                     if ( dstijk(1,j) /= PSMILe_Undef .and. &
!rr                          dstijk(2,j) /= PSMILe_Undef .and. &
!rr                          dstijk(1,j)-shape(1,1,ipart)+1 <= offset .and. &
!rr                          dstijk(3,j) /= PSMILe_Undef ) ncpl = ncpl + 1
                  enddo

                  nprev = nprev + len_cpl (ipart)

               endif

            enddo

            send_info%nloc = ncpl

            do i = 1, ndim_3d
               Allocate (tgt_cell(i)%vector(ncpl*tgt_corners(i)), stat = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = ncpl

                  ierror = PRISM_Error_Alloc
                  call psmile_error ( ierror, 'tgt_cell(i)%vector', &
                       ierrp, 2, __FILE__, __LINE__ )
                  return
               endif
            end do

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

         nprev = 0
         jj = 0

         do ipart = 1, search%npart

            if (len_cpl (ipart) > 0) then

               ibeg = nprev + 1
               iend = nprev + len_cpl (ipart)

               if (search%grid_type == PRISM_Irrlonlat_regvrt) then

                  len = shape(2,1,ipart)-shape(1,1,ipart)+1
                  do i = 1, ndim_2d

                     do jpart = ibeg, iend, 5000
!cdir vector loopcnt=5000
                     do j = jpart, min(iend,jpart+5000-1)
!rr                        if ( dstijk(1,j) /= PSMILe_Undef .and. dstijk(2,j) /= PSMILe_Undef ) then
                           n = (dstijk(2,j)-shape(1,2,ipart))*len + &
                                dstijk(1,j)-shape(1,1,ipart)+1
                           tgt_cell(i)%vector(0*ncpl+j) = coords(i,ipart)%vector(n)
                           tgt_cell(i)%vector(1*ncpl+j) = coords(i,ipart)%vector(n+1)
                           tgt_cell(i)%vector(2*ncpl+j) = coords(i,ipart)%vector(n+len+1)
                           tgt_cell(i)%vector(3*ncpl+j) = coords(i,ipart)%vector(n+len)
!rr                        else
!rr                           do nn = 0, nb_neighbors-1
!rr                              tgt_cell(i)%vector(nn*ncpl+j) = 0.0
!rr                           enddo
!rr                        endif
                     enddo
                     enddo

                  end do

                  do j = ibeg, iend
                     if ( dstijk(3,j) /= PSMILe_Undef ) then
                        n = dstijk(3,j)-shape(1,3,ipart)+1
                        do nn = 0, nb_neighbors-1
                           tgt_cell(3)%vector(nn*ncpl+j) = coords(3,ipart)%vector(n)
                        enddo
                     else
                        do nn = 0, nb_neighbors-1
                           tgt_cell(3)%vector(nn*ncpl+j) = 0.0
                        enddo
                     endif
                  enddo

                  if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
                     ierrp (1) = search%grid_type
                     ierror = PRISM_Error_Internal
                     call psmile_error ( ierror, 'unsupported interpolation', &
                          ierrp, 1, __FILE__, __LINE__ )
                     return     
                  endif

               else if ( search%grid_type == PRISM_Reglonlatvrt ) then

                  len = shape(2,1,ipart)-shape(1,1,ipart)+1

                  do jpart = ibeg, iend, 5000
!cdir vector loopcnt=5000
                  do j = jpart, min(iend,jpart+5000-1)

                     if ( dstijk(1,j) /= PSMILe_Undef .and. dstijk(2,j) /= PSMILe_Undef ) then
                        n1 = dstijk(1,j)-shape(1,1,ipart)
!rr                     n2 = dstijk(2,j)-shape(1,2,ipart)
                        n2 = (dstijk(2,j)-shape(1,2,ipart))*len
                        tgt_cell(1)%vector(0*ncpl+j) = coords(1,ipart)%vector(n1+1)
                        tgt_cell(1)%vector(1*ncpl+j) = coords(1,ipart)%vector(n1+2)
                        tgt_cell(1)%vector(2*ncpl+j) = coords(1,ipart)%vector(n1+2)
                        tgt_cell(1)%vector(3*ncpl+j) = coords(1,ipart)%vector(n1+1)
                        tgt_cell(2)%vector(0*ncpl+j) = coords(2,ipart)%vector(n2+1)
                        tgt_cell(2)%vector(1*ncpl+j) = coords(2,ipart)%vector(n2+1)
                        tgt_cell(2)%vector(2*ncpl+j) = coords(2,ipart)%vector(n2+len+1)
                        tgt_cell(2)%vector(3*ncpl+j) = coords(2,ipart)%vector(n2+len+1)
                      else
                        do nn = 0, nb_neighbors-1
                           tgt_cell(1)%vector(nn*ncpl+j) = 0.0
                           tgt_cell(2)%vector(nn*ncpl+j) = 0.0
                        enddo
                      endif
                  end do
                  end do

                  if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
                     do j = ibeg, iend
                        if ( dstijk(3,j) /= PSMILe_Undef ) then
                           n3 = dstijk(3,j)-shape(1,3,ipart)+1
                           do nn = 0, 3
                              tgt_cell(3)%vector(nn*ncpl+j) = coords(3,ipart)%vector(n3)
                           enddo
                           do nn = 4, nb_neighbors-1
                              tgt_cell(3)%vector(nn*ncpl+j) = coords(3,ipart)%vector(n3+1)
                           enddo
                        else
                           do nn = 0, nb_neighbors-1
                              tgt_cell(3)%vector(nn*ncpl+j) = 0.0
                           enddo
                        endif
                     end do
                  else
                     do j = ibeg, iend
                        if ( dstijk(3,j) /= PSMILe_Undef ) then
                           n3 = dstijk(3,j)-shape(1,3,ipart)+1
                           do nn = 0, nb_neighbors-1
                              tgt_cell(3)%vector(nn*ncpl+j) = coords(3,ipart)%vector(n3)
                           enddo
                        else
                           do nn = 0, nb_neighbors-1
                              tgt_cell(3)%vector(nn*ncpl+j) = 0.0
                           enddo
                        endif
                     enddo
                  endif

               else if (search%grid_type == PRISM_Gaussreduced_regvrt) then

               if (len_cpl (ipart) > 0) then

                  totsiz =  size(coords(1,ipart)%vector)
                  offset =  totsiz/2

                  len_cpl(ipart) = 0

                  do jpart = ibeg, iend, 5000
!cdir vector loopcnt=5000
                  do j = jpart, min(iend,jpart+5000-1)

                     if ( dstijk(1,j) /= PSMILe_Undef .and. &
                          dstijk(1,j)-shape(1,1,ipart)+1 <= offset ) then

!rr                       dstijk(2,j) /= PSMILe_Undef .and. &
!rr                       dstijk(1,j)-shape(1,1,ipart)+1 <= offset .and. &
!rr                       dstijk(3,j) /= PSMILe_Undef ) then

                        n = dstijk(1,j)-shape(1,1,ipart)+1
                        jj = jj + 1
                        len_cpl(ipart) = len_cpl(ipart) + 1
 
                        tgt_cell(1)%vector(0*ncpl+jj) = coords(1,ipart)%vector(n)
                        tgt_cell(1)%vector(1*ncpl+jj) = coords(1,ipart)%vector(n+offset)
                        tgt_cell(1)%vector(2*ncpl+jj) = coords(1,ipart)%vector(n+offset)
                        tgt_cell(1)%vector(3*ncpl+jj) = coords(1,ipart)%vector(n)
                        tgt_cell(2)%vector(0*ncpl+jj) = coords(2,ipart)%vector(n)
                        tgt_cell(2)%vector(1*ncpl+jj) = coords(2,ipart)%vector(n)
                        tgt_cell(2)%vector(2*ncpl+jj) = coords(2,ipart)%vector(n+offset)
                        tgt_cell(2)%vector(3*ncpl+jj) = coords(2,ipart)%vector(n+offset)

                        n3 = dstijk(3,j)-shape(1,3,ipart)+1
                        do nn = 0, nb_neighbors-1
                           tgt_cell(3)%vector(nn*ncpl+jj) = coords(3,ipart)%vector(n3)
                        enddo
!rr                     else
!rr                        ! For Gauss-reduced grids the msklocs vector
!rr                        ! runs over all iparts
!rr                        send_info%msklocs(1,1)%vector(j) = .false.
                     endif

                  enddo
                  enddo

                  if ( search%msg_intersections%requires_conserv_remap == PSMILe_Conserv3D ) then
                     ierrp (1) = search%grid_type
                     ierror = PRISM_Error_Internal
                     call psmile_error ( ierror, 'unsupported interpolation', &
                          ierrp, 1, __FILE__, __LINE__ )
                     return     
                  endif

               endif ! len_cpl > 0

               else

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

               endif ! search%grid_type 

               nprev = nprev + len_cpl (ipart)

            endif

         end do ! ipart

      else
!
!===> Interpolation using grid data ------------------------------------
!
         use_target_grid = .not. ( &
                        search%grid_type == PRISM_Irrlonlatvrt .and. &
                        search%npart == 1 .and. &
                        ncpl == (shape(2,1,1)-shape(1,1,1)+1)* &
                                (shape(2,2,1)-shape(1,2,1)+1)* &
                                (shape(2,3,1)-shape(1,3,1)+1) )

         if (use_target_grid) then
!
!===> ... Allocate vectors for PSMILe_Trs_set_tgt_epio3d_dble
!
         do i = 1, ndim_3d
            Allocate (tgt_grid(i)%vector(ncpl), stat = ierror)
            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = ncpl

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'tgt_grid(i)%vector', &
                    ierrp, 2, __FILE__, __LINE__ )
               return
            endif
         end do
!
!===> ... Extract coordinates of target (destination) points
!
         if ( search%grid_type == PRISM_Reglonlatvrt .or. &
              search%grid_type == PRISM_Irrlonlat_regvrt  .or. &
              search%grid_type == PRISM_Gaussreduced_regvrt) then
!
!           Search coordinates are a 2d array and 1d vector
!
            nprev = 0
            do ipart = 1, search%npart
               if (len_cpl (ipart) > 0) then
                  ibeg = nprev + 1
                  iend = nprev + len_cpl (ipart)
!
                  call psmile_extract_indices_2d_dble (                       &
                          coords(1,ipart)%vector, shape (:,1:ndim_2d, ipart), &
                          dstijk (:, ibeg:iend),  len_cpl(ipart),             &
                          tgt_grid(1)%vector(ibeg:iend), ierror)
                  if (ierror /= 0) return
!
                  call psmile_extract_indices_2d_dble (                       &
                          coords(2,ipart)%vector, shape (:,1:ndim_2d, ipart), &
                          dstijk (:, ibeg:iend),  len_cpl(ipart),             &
                          tgt_grid(2)%vector(ibeg:iend), ierror)
                  if (ierror /= 0) return
!
                  tgt_grid(3)%vector(ibeg:iend) = &
                     coords(3,ipart)%vector(dstijk(3,ibeg:iend)-shape(1,3,ipart)+1)

                  nprev = nprev + len_cpl (ipart)
               endif
            end do ! ipart
!
         else if (search%grid_type == PRISM_Irrlonlatvrt) then
!
!           Each search coordinate is stored in a 3-d array
!
            nprev = 0
            do ipart = 1, search%npart
               if (len_cpl (ipart) > 0) then
                  ibeg = nprev + 1
                  iend = nprev + len_cpl (ipart)
!
                  do i = 1, ndim_3d
                  call psmile_extract_indices_3d_dble (                  &
                              coords(i,ipart)%vector, shape(:,:, ipart), &
                              dstijk (:, ibeg:iend),  len_cpl(ipart),    &
                              tgt_grid(i)%vector(ibeg:iend),             &
                              ierror)
                  if (ierror > 0) return
                  end do

                  nprev = nprev + len_cpl (ipart)
               endif
            end do ! ipart
!
         else
!
!===> Something is missing !
!
            ierrp (1) = search%grid_type
            ierror = PRISM_Error_Internal
            call psmile_error ( ierror, 'unsupported grid type', &
                 ierrp, 1, __FILE__, __LINE__ )
            return          

         endif

         else
!
!===> ...... note: search%npart = 1
!
               do i = 1, ndim_3d
                  tgt_grid(i)%vector => coords(i,1)%vector
               end do
         endif

      endif ! use_target_cell, use_target_grid
!
!------------------------------------------------------------------------
!  Search indices of interpolation bases
!------------------------------------------------------------------------
!
!===> Allocate neighbour array (indices of interpolation bases)
!
      Allocate (neighbors_3d(ndim_3d, ncpl, nb_neighbors), STAT = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = ndim_3d * ncpl * nb_neighbors

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'neighbors_3d', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!===> Prepare data for extra and global search
!
      extra_search%global_search = global_search
!
      call psmile_neigh_extra_search_init (search, grid_id, extra_search, ierror)
      if (ierror > 0) return
!
!===> Perform search of indices
!
      search_mode = 3

      j = 1
         do while (j <= n_intmethods)
!
         select case (interpolation_methods(j))

!------------------------------------------------------------------------
!        Nearest neighbour in 3d
!------------------------------------------------------------------------

         case (PSMILe_nnghbr3D)
!
!===> ... Transform Longitudes and Latitudes coordinates
!         of source grid
!
         if (.not. transformed) then
            call psmile_info_coords_irreg2_dble (           &
                  mp%coords_pointer%coords_dble(1)%vector,  &
                  mp%coords_pointer%coords_dble(2)%vector,  &
                  mp%coords_pointer%coords_dble(3)%vector,  &
                  mp%coords_pointer%coords_shape,           &
                  grid_valid_shape,                         &
                  sinvec, cosvec, ierror)
            if (ierror > 0) return

            transformed = .true.
         endif
!
!====> ... Look for nearest neighbours
!
         nprev = 0

         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then
            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_near_irr2_3d_dble (grid_id,  &
                  tgt_grid(1)%vector,                         &
                  tgt_grid(2)%vector,                         &
                  tgt_grid(3)%vector,                         &
                  mp%coords_pointer%coords_dble(1)%vector,    &
                  mp%coords_pointer%coords_dble(2)%vector,    &
                  mp%coords_pointer%coords_dble(3)%vector,    &
                  mp%coords_pointer%coords_shape,             &
                  mask_array, mask_shape, use_mask,           &
                  sinvec%vector, cosvec%vector,               &
                  grid_valid_shape, search_mode,              &
                  send_info%srclocs(1,1)%vector,              &
                  ncpl, nprev, len_cpl (ipart),               &
                  neighbors_3d, nb_neighbors,                 &
                  extra_search, ierror)
                  if (ierror > 0) return

                  nprev = nprev + len_cpl (ipart)
            end if
            end do ! ipart

         else
            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_near_irreg2_dble (grid_id,  &
                  tgt_grid(1)%vector(nprev+1:),              &
                  tgt_grid(2)%vector(nprev+1:),              &
                  tgt_grid(3)%vector(nprev+1:),              &
                  mp%coords_pointer%coords_dble(1)%vector,   &
                  mp%coords_pointer%coords_dble(2)%vector,   &
                  mp%coords_pointer%coords_dble(3)%vector,   &
                  mp%coords_pointer%coords_shape,            &
                  mask_array, mask_shape, use_mask,          &
                  sinvec%vector, cosvec%vector,              &
                  grid_valid_shape, search_mode,             &
                  send_info%srclocs(indl, ipart)%vector,     &
                  send_info%srclocs(indz, ipart)%vector,     &
                  send_info%npoints(1:indz, ipart),          &
                  ncpl, nprev,                               &
                  neighbors_3d, nb_neighbors,                &
                  extra_search, ierror)

               if (ierror > 0) return

               nprev = nprev + len_cpl (ipart)
            end if 
            end do ! ipart
         endif
!
         case (PSMILe_nnghbr2D)
!
!===> ... Transform Longitudes and Latitudes coordinates of source grid
!
         if (.not. transformed) then
            call psmile_info_coords_irreg2_dble (          &
                  mp%coords_pointer%coords_dble(1)%vector, &
                  mp%coords_pointer%coords_dble(2)%vector, &
                  mp%coords_pointer%coords_dble(3)%vector, &
                  mp%coords_pointer%coords_shape,          &
                  grid_valid_shape,                        &
                  sinvec, cosvec, ierror)
               if (ierror > 0) return

            transformed = .true.
         endif
!
!====> ... Look for nearest neighbours
!
         nprev = 0
         search_mode = 2

         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then
            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_near_irr2_3d_dble (grid_id,  &
                  tgt_grid(1)%vector,                         &
                  tgt_grid(2)%vector,                         &
                  tgt_grid(3)%vector,                         &
                  mp%coords_pointer%coords_dble(1)%vector,    &
                  mp%coords_pointer%coords_dble(2)%vector,    &
                  mp%coords_pointer%coords_dble(3)%vector,    &
                  mp%coords_pointer%coords_shape,             &
                  mask_array, mask_shape, use_mask,           &
                  sinvec%vector, cosvec%vector,               &
                  grid_valid_shape, search_mode,              &
                  send_info%srclocs(1,1)%vector,              &
                  ncpl, nprev, len_cpl (ipart),               &
                  neighbors_3d, nb_neighbors,                 &
                  extra_search, ierror)
                  if (ierror > 0) return

                  nprev = nprev + len_cpl (ipart)
            end if
            end do ! ipart

         else
            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_near_irreg2_dble (grid_id,  &
                  tgt_grid(1)%vector(nprev+1:),              &
                  tgt_grid(2)%vector(nprev+1:),              &
                  tgt_grid(3)%vector(nprev+1:),              &
                  mp%coords_pointer%coords_dble(1)%vector,   &
                  mp%coords_pointer%coords_dble(2)%vector,   &
                  mp%coords_pointer%coords_dble(3)%vector,   &
                  mp%coords_pointer%coords_shape,            &
                  mask_array, mask_shape, use_mask,          &
                  sinvec%vector, cosvec%vector,              &
                  grid_valid_shape, search_mode,             &
                  send_info%srclocs(indl, ipart)%vector,     &
                  send_info%srclocs(indz, ipart)%vector,     &
                  send_info%npoints(1:indz, ipart),          &
                  ncpl, nprev,                               &
                  neighbors_3d, nb_neighbors,                &
                  extra_search, ierror)

               if (ierror > 0) return

               nprev = nprev + len_cpl (ipart)
            end if
            end do ! ipart
         endif

!-----------------------------------------------------------------------
!     Linear, Bi-linear, Bi-linear/Linear or Tri-linear interpolation
!-----------------------------------------------------------------------

      case (PSMILe_trilinear, PSMILe_bilinear, PSMILe_bicubic) 

!        interpolation in 3d
         interpolation_mode = 3
         search_mode = 3
!
         if (interpolation_methods(j) == PSMILe_bilinear .or. &
             interpolation_methods(j) == PSMILe_bicubic) then
            if (interpolation_methods(j+1) == PSMILe_linear) then
               j = j + 1
!
!   Should not this be search_mode = 2 and we have to find 
!   the 4 nearest neighbours two times?
!
               search_mode = 3
            else
!              only bi-linear (2d) interpolation
               interpolation_mode = 2
               search_mode = 2
            endif
         endif
!
         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then

            if ( cubic ) then
               call psmile_neigh_tricu_3d (                             &
                    grid_valid_shape,                                   &
                    interpolation_mode, grid_info%cyclic,               &
                    send_info%srclocs(1,1)%vector,                      &
                    ncpl, neighbors_3d, nb_neighbors,                   &
                    ierror)
            else
               call psmile_neigh_trili_3d (                             &
                    grid_valid_shape,                                   &
                    interpolation_mode, grid_info%cyclic,               &
                    send_info%srclocs(1,1)%vector,                      &
                    ncpl, neighbors_3d, nb_neighbors,                   &
                    ierror)
            endif
            if (ierror /= 0) return

         else
            nprev = 0

            do ipart = 1, search%npart
               if (len_cpl (ipart) > 0) then

                  if ( cubic ) then
                     call psmile_neigh_tricu_irreg2 (                    &
                          grid_valid_shape,                              &
                          interpolation_mode, grid_info%cyclic,          &
                          send_info%srclocs(indl, ipart)%vector,         &
                          send_info%srclocs(indz, ipart)%vector,         &
                          send_info%npoints(1:indz, ipart),              &
                          ncpl, nprev,                                   &
                          neighbors_3d, nb_neighbors, ierror)
                  else
                     call psmile_neigh_trili_irreg2 (                    &
                          grid_valid_shape,                              &
                          interpolation_mode, grid_info%cyclic,          &
                          send_info%srclocs(indl, ipart)%vector,         &
                          send_info%srclocs(indz, ipart)%vector,         &
                          send_info%npoints(1:indz, ipart),              &
                          ncpl, nprev,                                   &
                          neighbors_3d, nb_neighbors, ierror)
                  endif

                  if (ierror /= 0) return

                  nprev = nprev + len_cpl (ipart)
               end if
            end do ! ipart

         endif

!-----------------------------------------------------------------------
!     Conservative Remapping
!-----------------------------------------------------------------------
!
! needs to be redesigned.
!
      case (PSMILe_conserv2D)

         Allocate(nbr_cells(ncpl))

!        interpolation in 3d
         interpolation_mode = 3
!
         corner_shape_3d = 1

         do i = 1, ndim_2d
            corner_shape_3d(:,i,1) = cp%corner_shape(:,i)
            corner_shape_3d(:,i,2) = cp%corner_shape(:,i)
         enddo

         corner_shape_3d(:,i,3) = cp%corner_shape(:,3)

         if (interpolation_methods(j+1) == PSMILe_linear) then
            j = j + 1
         else
!           only 2d conservative interpolation
            interpolation_mode = 2
            search_mode = 2
         endif

         ipart = search%npart

         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then

            nbr_corners(1:2) = 4
            nbr_corners(3)   = 2

            call psmile_neigh_cells_3d_dble ( use_how,      &
               grid_valid_shape,                            &
               interpolation_mode,  grid_info%cyclic,       &
               corner_shape_3d, nbr_corners,                &
               cp%corners_dble(1)%vector,                   &
               cp%corners_dble(2)%vector,                   &
               cp%corners_dble(3)%vector,                   &
               search, control, tgt_cell, tgt_corners,      &
               send_info%npoints(1,1),                      &
               send_info%srclocs(1,1)%vector,               &
               send_info%msklocs(1,1)%vector,               &
               ncpl, nb_neighbors, nbr_cells, ierror )
            if (ierror /= 0) return

         else if (send_info%nvec == 2) then
!        else if (search%grid_type == PRISM_Irrlonlat_Regvrt .or. &
!                 search%grid_type == PRISM_Gaussreduced_Regvrt) then

            nbr_corners(1:2) = 4
            nbr_corners(3)   = 2

            call psmile_neigh_cells_irreg2_dble ( use_how,   &
                 grid_valid_shape, interpolation_mode,       &
                 grid_info%cyclic,                           &
                 corner_shape_3d, nbr_corners,               &
                 cp%corners_dble(1)%vector,                  &
                 cp%corners_dble(2)%vector,                  &
                 search, control, tgt_cell, tgt_corners(1),  &
                 send_info%npoints(1:indz, 1:ipart),         &
                 send_info%srclocs(1:indz, 1:ipart),         &
                 send_info%msklocs(1:indz, 1:ipart),         &
                 ncpl, nb_neighbors, nbr_cells, ierror )

            if (ierror /= 0) return

         else if (send_info%nvec == ndim_3d) then

            if (send_info%nloc > 0) then
               call psmile_neigh_cells_3d_reg_dble (       &
                    grid_valid_shape,                      &
                    interpolation_mode, grid_info%cyclic,  &
                    grid_id, search, coords,               &
                    send_info%npoints(1:ndim_3d, 1:ipart), &
                    send_info%srclocs(1:ndim_3d, 1:ipart), &
                    ncpl, nbr_cells, ierror )

               if (ierror /= 0) return
            end if

         endif

         if ( associated(search%boundary_cell) .and. global_cell_search ) then

            nprev = 0
            nn = 1

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

            do ipart = 1, search%npart

               if (len_cpl (ipart) > 0) then

                  ibeg = nprev + 1
                  iend = nprev + len_cpl (ipart)

                  if ( search%grid_type == PRISM_Irrlonlat_regvrt ) then

                     leni = shape(2,1,ipart)-shape(1,1,ipart)+1

                     ! We do not have to loop over the whole j. Instead
                     ! we can loop over the size of boundary_cell(:,1)
                     ! and take j=boundary_cell(:,1). In this case
                     ! we have to determine the correct ipart for each j.

                     do jpart = ibeg, iend, 5000
                        do j = jpart, min(iend,jpart+5000-1)

                           if ( dstijk(1,j) /= PSMILe_Undef .and. &
                                dstijk(2,j) /= PSMILe_Undef .and. &
                                nn <= len ) then

                              n = (dstijk(2,j)-shape(1,2,ipart))*leni + &
                                   dstijk(1,j)-shape(1,1,ipart)+1

                              if ( j == search%boundary_cell(nn,1) ) then
                                 n2 = ( n-1) / leni + 1
                                 n1 = n - (n2-1)*leni
                                 search%boundary_cell(nn,2) = &
                                     n1 + search%global_index(ipart)%vector(1)
                                 search%boundary_cell(nn,3) = &
                                     n2 + search%global_index(ipart)%vector(3)
                                 search%boundary_cell(nn,4) = &
                                     dstijk(3,j)-shape(1,3,ipart)+1 + &
                                          search%global_index(ipart)%vector(5)
                                 nn = nn + 1
                              endif

                           endif
                        enddo
                     enddo

                  else if ( search%grid_type == PRISM_Reglonlatvrt ) then

                     ! We do not have to loop over the whole j. Instead
                     ! we can loop over the size of boundary_cell(:,1)
                     ! and take j=boundary_cell(:,1). In this case
                     ! we have to determine the correct ipart for each j.

                     do jpart = ibeg, iend, 5000
                        do j = jpart, min(iend,jpart+5000-1)

                           if ( dstijk(1,j) /= PSMILe_Undef .and. &
                                dstijk(2,j) /= PSMILe_Undef .and. &
                                nn <= len ) then

                              if ( j == search%boundary_cell(nn,1) ) then

                                 search%boundary_cell(nn,2) = dstijk(1,j)-shape(1,1,ipart)+1 + &
                                                              search%global_index(ipart)%vector(1)
                                 search%boundary_cell(nn,3) = dstijk(2,j)-shape(1,2,ipart)+1 + &
                                                              search%global_index(ipart)%vector(3)
                                 search%boundary_cell(nn,4) = dstijk(3,j)-shape(1,3,ipart)+1 + &
                                                              search%global_index(ipart)%vector(5)
                                 nn = nn + 1
                              endif

                           endif
                        enddo
                     enddo

                  else if ( search%grid_type == PRISM_Gaussreduced_regvrt ) then

                     ! We do not have to loop over the whole j. Instead
                     ! we can loop over the size of boundary_cell(:,1)
                     ! and take j=boundary_cell(:,1). In this case
                     ! we have to determine the correct ipart for each j.

                     do jpart = ibeg, iend, 5000
                        do j = jpart, min(iend,jpart+5000-1)

                           if ( dstijk(1,j) /= PSMILe_Undef .and. &
                                dstijk(2,j) /= PSMILe_Undef .and. &
                                nn <= len ) then

                              if ( j == search%boundary_cell(nn,1) ) then

                                 search%boundary_cell(nn,2) = dstijk(1,j)-shape(1,1,ipart)+1 + &
                                                              search%global_index(ipart)%vector(1)
                                 search%boundary_cell(nn,3) = 1
                                 search%boundary_cell(nn,4) = dstijk(3,j)-shape(1,3,ipart)+1 + &
                                                              search%global_index(ipart)%vector(5)
                                 nn = nn + 1
                              endif

                           endif
                        enddo
                     enddo

                  else
                     
                     !
                     !===> Something is missing !
                     !
                     ierrp (1) = search%grid_type
                     ierror = PRISM_Error_Internal
                     call psmile_error ( ierror, 'unsupported grid type', &
                          ierrp, 1, __FILE__, __LINE__ )
                     return

                  endif

                  nprev = iend

               endif

            enddo

         endif

         if ( global_cell_search ) &
            call psmile_global_search_cell_dble ( grid_id, var_id,  &
                  comp_info, send_info, search, extra_search, ncpl, &
                  nbr_cells, n_intmethods, interpolation_methods,   &
                  interpolation_search, ierror )

         if ( associated(search%boundary_cell) .and. global_cell_search ) then
            !
            ! Erase information about target cells which have been
            ! transferred to other source processes and will be handled from there.
            !

            new_len_cpl = len_cpl
            new_ncpl    = send_info%nloc

            if ( new_ncpl < ncpl ) then
               !
               ! Allocate memory for reduced target information
               !
               allocate(new_nbr_cells(new_ncpl), new_dstijk(3,new_ncpl), stat = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = send_info%nloc

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

               do i = 1, ndim_3d
                  Allocate (new_tgt_cell(i)%vector(new_ncpl*tgt_corners(i)), stat = ierror)
                  if ( ierror > 0 ) then
                     ierrp (1) = ierror
                     ierrp (2) = send_info%nloc*tgt_corners(i)

                     ierror = PRISM_Error_Alloc
                     call psmile_error ( ierror, 'new_tgt_cell(i)%vector', &
                          ierrp, 2, __FILE__, __LINE__ )
                     return
                  endif
               end do

               nn = 1

               ipart = 1

               ibeg = 1
               iend = len_cpl (ipart)

               nprev = iend
               !
               ! Extract target information and store in reduced (new) list
               !
               do n = 1, ncpl

                  if ( n > iend ) then
                     ipart = ipart + 1
                     ibeg = nprev + 1
                     iend = nprev + len_cpl (ipart)
                     nprev = iend
                  endif
!rr                  do while ( n > iend )
!rr                     ipart = ipart + 1
!rr                     ibeg = nprev + 1
!rr                     iend = nprev + len_cpl (ipart)
!rr                     nprev = iend
!rr                  enddo
                  !
                  ! nbr_cells(n) have been set to zero in psmile_gobal_search_cell
                  ! for target cells treated elsewhere.
                  !
                  if ( nbr_cells(n) > 0 ) then
                     new_nbr_cells(nn) = nbr_cells(n)
                     new_dstijk(:,nn)  = dstijk(:,n)
                     nn = nn + 1
                  else
                     new_len_cpl(ipart) = new_len_cpl(ipart) - 1
                  endif

               enddo

               len_cpl = new_len_cpl

               do i = 1, ndim_2d
               nn = 1
               do n = 1, ncpl
                  if ( nbr_cells(n) > 0 ) then
#ifdef DEBUGX
                     print *, ' - shift tgt cell idx ', n , ' -> ', nn
#endif
                     new_tgt_cell(i)%vector(0*new_ncpl+nn) = tgt_cell(i)%vector(0*ncpl+n)
                     new_tgt_cell(i)%vector(1*new_ncpl+nn) = tgt_cell(i)%vector(1*ncpl+n)
                     new_tgt_cell(i)%vector(2*new_ncpl+nn) = tgt_cell(i)%vector(2*ncpl+n)
                     new_tgt_cell(i)%vector(3*new_ncpl+nn) = tgt_cell(i)%vector(3*ncpl+n)
                     nn = nn + 1
                  endif
               enddo
               enddo

               do n = 1, ncpl
                  nn = 1
                  if ( nbr_cells(n) > 0 ) then
                     do n3 = 0, nb_neighbors-1
                        new_tgt_cell(3)%vector(n3*new_ncpl+nn) = tgt_cell(3)%vector(n3*ncpl+n)
                     enddo
                     nn = nn + 1
                  endif
               enddo
               !
               ! Deallocate old memory
               !
               deallocate(nbr_cells, send_info%dstijk)
               do i = 1, ndim_3d
                 deallocate(tgt_cell(i)%vector)
               enddo

               Allocate(send_info%dstijk(3,new_ncpl), stat = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = send_info%nloc

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

               nbr_cells          => new_nbr_cells

               tgt_cell(1)%vector => new_tgt_cell(1)%vector
               tgt_cell(2)%vector => new_tgt_cell(2)%vector
               tgt_cell(3)%vector => new_tgt_cell(3)%vector

               send_info%dstijk = new_dstijk

               Deallocate(new_dstijk)

               ncpl = send_info%nloc

            endif

         endif

#ifdef TODO
      case (PSMILe_conserv3D)

!        interpolation in 3d
         interpolation_mode = 3
         search_mode = 3

         corner_shape_3d = 1

         if (search%grid_type == PRISM_Irrlonlat_regvrt) then
            do i = 1, ndim_2d
               corner_shape_3d(:,i,1) = cp%corner_shape(:,i)
               corner_shape_3d(:,i,2) = cp%corner_shape(:,i)
            enddo
            corner_shape_3d(:,i,3) = cp%corner_shape(:,3)
         else if (search%grid_type == PRISM_Irrlonlatvrt) then
            do i = 1, ndim_3d
               corner_shape_3d(:,:,i) = cp%corner_shape(:,:)
            enddo
         else if (search%grid_type == PRISM_Reglonlatvrt) then
            do i = 1, ndim_3d
               corner_shape_3d(:,i,i) = cp%corner_shape(:,i)
            enddo
         else if (search%grid_type == PRISM_Gaussreduced_regvrt) then
            do i = 1, ndim_3d
               corner_shape_3d(:,i,i) = cp%corner_shape(:,i)
            enddo
         else
!
!===> Something is missing !
!
            ierrp (1) = search%grid_type
            ierror = PRISM_Error_Internal
            call psmile_error ( ierror, 'unsupported grid type', &
                 ierrp, 1, __FILE__, __LINE__ )
            return

         endif

         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then
            
            nbr_corners(1:2) = 4
            nbr_corners(3)   = 2

            call psmile_neigh_cells_3d_dble ( use_how,  &
               grid_valid_shape,                        &
               interpolation_mode,  grid_info%cyclic,   &
               corner_shape_3d, nbr_corners,            &
               cp%corners_dble(1)%vector,               &
               cp%corners_dble(2)%vector,               &
               cp%corners_dble(3)%vector,               &
               search, control, tgt_cell, tgt_corners,  &
               send_info%npoints(1,1),                  &
               send_info%srclocs(1,1)%vector,           &
               send_info%msklocs(1,1)%vector,           &
               ncpl, nb_neighbors, nbr_cells, ierror )
            if (ierror /= 0) return

         else if (send_info%nvec == 2) then
!        else if (search%grid_type == PRISM_Irrlonlat_Regvrt .or. &
!                 search%grid_type == PRISM_Gaussreduced_Regvrt) then

            nbr_corners = 4

            call psmile_neigh_cells_irreg2_dble ( use_how,  &
                 grid_valid_shape, interpolation_mode,      &
                 grid_info%cyclic,                          &
                 corner_shape_3d, nbr_corners,              &
                 cp%corners_dble(1)%vector,                 &
                 cp%corners_dble(2)%vector,                 &
                 search, control, tgt_cell, tgt_corners(1), &
                 send_info%npoints(1:indz, 1:ipart),        &
                 send_info%srclocs(1:indz, 1:ipart),        &
                 send_info%msklocs(1:indz, 1:ipart),        &
                 ncpl, nb_neighbors, nbr_cells, ierror )

            if (ierror /= 0) return

         else ! if (send_info%nvec == ndim_3d) then

            if (send_info%nloc > 0) then
               call psmile_neigh_cells_3d_reg_dble (       &
                    grid_valid_shape,                      &
                    interpolation_mode, grid_info%cyclic,  &
                    grid_id, search, coords,               &
                    send_info%npoints(1:ndim_3d, 1:ipart), &
                    send_info%srclocs(1:ndim_3d, 1:ipart), &
                    ncpl, nbr_cells, ierror )

               if (ierror /= 0) return
            end if

         endif
#endif

#ifdef TODO
      case (PSMILe_user3D)
         call psmile_trs_give_
         if (ierror /= 0) return
#endif

!     Error: unsupported interpolation method

      case default
         ierrp (1) = interpolation_methods (j)
         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'unsupported 3d interpolation method', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      end select
!
!===> ... Next interpolation
!
      j = j + 1
      end do ! while ( j < n_intmethods )
!
#ifdef DEBUG_TRACE
      if ( use_target_grid ) then
!
!cdir vector
         do ictl = 1, ncpl
         if (dstijk (1, ictl) == ictl_ind (1) .and. &
             dstijk (2, ictl) == ictl_ind (2) .and. &
             dstijk (3, ictl) == ictl_ind (3)) exit
         end do
!
      if (ictl <= ncpl) then
         print *, 'vor global search: ictl, dstijk (:, i)', &
                  ictl, ';', dstijk (:, ictl)
         print *, 'ictl target coord', ictl, tgt_grid(1)%vector(ictl), &
                                             tgt_grid(2)%vector(ictl), &
                                             tgt_grid(3)%vector(ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      else
         print *, 'psmile_info_trs_irreg2_dble: ictl_ind not found', ictl_ind
      endif
      endif
#endif
!
!
!===> Perform global search if required
!
      if (global_search) then
         if (interpolation_methods(1) /= PSMILe_nnghbr3D .and. &
             interpolation_methods(1) /= PSMILe_nnghbr2D) then
!
!           ... Global Search for non-nearest neighbour
!
            call psmile_neigh_global_points (search, extra_search,        &
                  mask_array, mask_shape, src_mask_available, use_how(1), &
                  grid_id, grid_valid_shape,                           &
                  neighbors_3d, ncpl, nb_neighbors, len_cpl, ierror)
            if (ierror > 0) return

            if (extra_search%n_req > 0) then
!
               call psmile_global_search_dble (comp_info,                 &
                         control, len_cpl,                                &
                         var_id, grid_valid_shape, search, tgt_grid,      &
                         neighbors_3d, ncpl, nb_neighbors, extra_search,  &
                         interpolation_methods, interpolation_search,     &
                         n_intmethods,                                    &
                         send_index, src_mask_available, use_mask,        &
                         use_how(1), PRISM_Irrlonlat_Regvrt, ierror)
               if (ierror > 0) return
            endif
!
         else
!
!           ... Global Search for nearest neighbour
!
         endif
      endif
!
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'after global_search: i, dstijk (:, i)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
!
!===> Get points of extra nearest neighbour search and
!     mask out points.
!
      if (interpolation_methods(1) /= PSMILe_nnghbr3D) then
         call psmile_neigh_extra_points (search, extra_search,           &
                 mask_array, mask_shape, src_mask_available, use_how,    &
                 grid_valid_shape,                                       &
                 neighbors_3d, ncpl, nb_neighbors, len_cpl, ierror)
         if (ierror > 0) return
      endif
!
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'after neigh_extra: ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
!
!===> Perform nearest neighbour search for extra points which were not found 
!
      if (extra_search%n_extra > 0) then
!
!===> ... Allocate distances for global search
!         in order to store distances of local nearest neighbour search.
!
         if (global_search) then
            call psmile_neigh_extra_search_dble (search, extra_search,  &
                    nb_extra, ierror)
            if (ierror > 0) return
         endif
!
!        Search mode for extra search is set to 2d search for 
!        irrlonlat_regvrt target grids 
!
         search_mode = 2
!
         nprev = 0
!
!===> ... Transform Longitudes and Latitudes coordinates if necessary
!
         if (.not. transformed) then
            call psmile_info_coords_irreg2_dble (                       &
                  mp%coords_pointer%coords_dble(1)%vector,              &
                  mp%coords_pointer%coords_dble(2)%vector,              &
                  mp%coords_pointer%coords_dble(3)%vector,              &
                  mp%coords_pointer%coords_shape,                       &
                  grid_valid_shape,                                     &
                  sinvec, cosvec, ierror)
            if (ierror > 0) return

            transformed = .true.
         endif
!
!===> ... Search
!
         if (send_info%nvec == 1) then
!        if (search%grid_type == PRISM_Irrlonlatvrt) then
            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_nearx_irr2_3d_dble (grid_id,     &
                     tgt_grid(1)%vector,                          &
                     tgt_grid(2)%vector,                          &
                     tgt_grid(3)%vector,                          &
                     mp%coords_pointer%coords_dble(1)%vector,     &
                     mp%coords_pointer%coords_dble(2)%vector,     &
                     mp%coords_pointer%coords_dble(3)%vector,     &
                     mp%coords_pointer%coords_shape,              &
                     mask_array, mask_shape, src_mask_available,  &
                     sinvec%vector, cosvec%vector,                &
                     grid_valid_shape, search_mode,               &
                     send_info%srclocs(1,1)%vector,               &
                     ncpl, nprev, len_cpl (ipart),                &
                     neighbors_3d, nb_extra,                      &
                     extra_search, ierror)

               if (ierror > 0) return

               nprev = nprev + len_cpl (ipart)
            end if
            end do ! ipart
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'after PSMILe_Neigh_nearx_irr2_3d_dble'
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
         else

            do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then

               call psmile_neigh_nearx_irreg2_dble (grid_id,      &
                     tgt_grid(1)%vector, tgt_grid(2)%vector,      &
                     tgt_grid(3)%vector,                          &
                     mp%coords_pointer%coords_dble(1)%vector,     &
                     mp%coords_pointer%coords_dble(2)%vector,     &
                     mp%coords_pointer%coords_dble(3)%vector,     &
                     mp%coords_pointer%coords_shape,              &
                     mask_array, mask_shape, src_mask_available,  &
                     sinvec%vector, cosvec%vector,                &
                     grid_valid_shape, search_mode,               &
                     send_info%srclocs(indl, ipart)%vector,       &
                     send_info%srclocs(indz, ipart)%vector,       &
                     send_info%npoints(1:indz, ipart),            &
                     ncpl, nprev,                                 &
                     neighbors_3d, nb_extra,                      &
                     extra_search, ierror)

               if (ierror > 0) return

               nprev = nprev + len_cpl (ipart)
            end if
            end do ! ipart
         endif
!
!
!===> Perform global search of nearest neighbour points
!
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'before PSMILe_global_search_nnx_dble'
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
         if (global_search) then
            call psmile_global_search_nnx_dble (comp_info,    &
                  search, var_id,                             &
                  tgt_grid(1)%vector,                         &
                  tgt_grid(2)%vector,                         &
                  tgt_grid(3)%vector,                         &
                  neighbors_3d, ncpl, nb_neighbors, nb_extra, &
                  extra_search, send_index, ierror)
            if (ierror > 0) return
         endif
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'after PSMILe_global_search_nnx_dble'
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
!
!===> ... Deallocate indices to be searched within the extra sweep
!
         call psmile_neigh_extra_search_clean (search, extra_search, ierror)
      endif
!
!===> Deallocate vectors for termporary values
!     Should we save for another search on this method?
!
      if (transformed) then
         Deallocate (cosvec%vector)
         Deallocate (sinvec%vector)
      endif
!
!-----------------------------------------------------------------------------
!  Transform 3d indices into 1d locations and
!  create vectors required by PSMILe_Trs_set_src_epio3d_dble.
!-----------------------------------------------------------------------------

      if ( use_source_cell ) then

         nbr_cells_tot = sum(nbr_cells) ! already contains num2recv
         nb_corners = nb_neighbors ! 4
!
!===> Check for empty lists as the neighcells_3d may have been driven
!      down to zero due to the transfer of cells to neighboring processes.
!
         if ( nbr_cells_tot == 0 ) then
            send_index = PRISM_UNDEFINED
            go to 1100
         endif

         Allocate (neighcells(nbr_cells_tot, nb_corners), &
                   source_cell_index(nbr_cells_tot), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = nbr_cells_tot * ( nb_corners + 1 )
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'neighcells & source_cell_index', &
                 ierrp, 2, __FILE__, __LINE__ )
            return
         endif

!rr         corner_shape(1,:) = grid_valid_shape(1,:)
!rr         corner_shape(2,1) = grid_valid_shape(2,1)
!rr         corner_shape(2,2) = grid_valid_shape(2,2)
!rr         corner_shape(2,3) = grid_valid_shape(2,3)
         !
         ! only work on reference index for cell
         !
         nb_corners = 1

         call psmile_compact_neighbors_3d (neighcells_3d, nbr_cells_tot,  &
                                           nb_corners, grid_valid_shape,  &
                                           extra_search, send_info,       &
                                           neighcells, ierror)
#ifdef move0_for_cells 
!
!===> ... Move 0 (dummy) entries to end of list
!         This is probably not necessary as for each cell the indices
!         do not contain any dummy entries. This could be used however
!         to move complete cells out. Don't do this. neighcells may contain
!         zeros due to global marker for global search.
!
         call psmile_move0_neighbors (neighcells, nbr_cells_tot, nb_corners, &
                ierror)
         if (ierror > 0) return
#endif /* move0_for_cells */
!
!===>   Build source cell array using n_list and list_entries
!
         src_cell_size = send_info%n_list * nb_neighbors

         do i = 1, ndim_3d
            Allocate (src_cell(i)%vector(src_cell_size), stat = ierror)
            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = src_cell_size

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'src_cell(i)%vector', &
                    ierrp, 2, __FILE__, __LINE__ )
               return
            endif
         end do
!
!===> Compact source corners and create new list_entries and n_list which is
!       used for the mask and later for the data
!
         corner_shape(1:2,1:2) = cp%corner_shape(1:2,1:2)
         corner_shape(:,3)     = mp%coords_pointer%coords_shape(:,3)

         call psmile_ccompact_irreg2_dble (                      &
                        send_info, grid_valid_shape,             &
                        corner_shape, cp%nbr_corners/2,          &
                        cp%corners_dble(1)%vector,               &
                        cp%corners_dble(2)%vector,               &
                        mp%coords_pointer%coords_dble(3)%vector, &
                        extra_search,                            &
                        src_cell_size, nbr_cells_tot,            &
                        source_cell_index,                       &
                        neighcells,                              &
                        src_cell(1)%vector,                      &
                        src_cell(2)%vector,                      &
                        src_cell(3)%vector, ierror )
         if (ierror /= 0) return

         if (associated(neighcells_3d)) Deallocate (neighcells_3d)

      else ! we have to work on grid points
!
!===> Allocate neighbour array (indices of interpolation bases)
!
         Allocate (neighbors(ncpl, nb_neighbors), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = ncpl * nb_neighbors

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'neighbors', &
                 ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'before compact neighbors_3d'
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
         call psmile_compact_neighbors_3d (neighbors_3d, ncpl,             &
                                           nb_neighbors, grid_valid_shape, &
                                           extra_search, send_info,        &
                                           neighbors, ierror)
         if (ierror /= 0) return
#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'after compact neighbors_3d'
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
      endif
#endif
!
!===> ... Move 0 (dummy) entries to end of list
!
         if (interpolation_methods(1) /= PSMILe_nnghbr3D) then
            call psmile_move0_neighbors (neighbors, ncpl, nb_neighbors, &
                                         ierror)
            if (ierror > 0) return
         endif

#ifdef DEBUG_TRACE
      if (ictl <= ncpl .and. use_target_grid) then
         print *, 'ictl, dstijk (:, ictl)', ictl, dstijk (:, ictl)
!
         do i = 1, nb_neighbors
            print *, 'neighbors_3d (:, ictl)', neighbors_3d (:, ictl, i)
         end do
!
         print *, 'neighbors (ictl)', neighbors (ictl, 1:nb_neighbors)
         if (allocated (ictl_neighbours)) deallocate (ictl_neighbours)
         allocate (ictl_neighbours(nb_neighbors))
         ictl_neighbours (1:nb_neighbors) = neighbors (ictl, 1:nb_neighbors)
      endif
#endif
         Deallocate (neighbors_3d, STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = ndim_3d*ncpl*nb_neighbors

            ierror = PRISM_Error_Dealloc
            call psmile_error ( ierror, 'neighbors_3d', &
                 ierrp, 2, __FILE__, __LINE__ )
            return
         endif

      endif ! use_source_cell
!
!=============================================================================
!  Transfer data to the coupler process
!=============================================================================
!
!  global comp id in dest process =
!        all_comp_infos(search%msg_intersections%all_comp_infos_comp_idx)%global_comp_id
!
!----------------------------------------------------------------------------
!     Transfer data on source 
!     Note: "list_entries" is generated by PSMILe_compact_neighbors_3d
!           if send_info%send_entire_valid_shape is not .true.
!----------------------------------------------------------------------------
!
!  id_src_mask > 0: src_mask is transferred (only for conservative)
!              = 0: src_mask is not transferred
!
      if ( use_target_cell ) then
         if ( send_info%nvec == 1 ) then
              send_info%npoints(1,1) = sum(len_cpl(1:search%npart))
         else if ( send_info%nvec == 2 ) then
            do ipart = 1, search%npart
               if ( send_info%npoints(indz,ipart) /= 0 ) then
               send_info%npoints(indl,ipart) = len_cpl (ipart)/ &
               send_info%npoints(indz,ipart)
               else
                  send_info%npoints(indl,ipart) = 0
               endif
            enddo
         endif
      endif

      src_size = send_info%n_list
      local_src_size = src_size - send_info%num2recv
!
!===> Check for empty lists
!
      if ( src_size == 0 ) then
         send_index = PRISM_UNDEFINED

         go to 1000
      endif
!
!===> For conservative remapping we transfer the source mask to the transformer
!
!rr      do j = 1, n_intmethods
!rr         if ( interpolation_methods(j) == PSMILe_conserv2D .or. &
!rr              interpolation_methods(j) == PSMILe_conserv3D ) then
!rr            if ( src_mask_available ) then
!rr               id_src_mask = 1
!rr               use_mask    = .false.
!rr            endif
!rr         endif
!rr      enddo
!
!===> Allocate and extract mask values if required
!
      if ( use_source_cell .and. src_mask_available) then
         use_mask = .false.
         id_src_mask = 1
!
         Allocate (src_mask(src_size), stat = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = src_size

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'src_mask', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif 
!
!===> ... Extract mask values of compact list
!         !!! mask of transformer should be also LOGICAL !!!
!
         call psmile_ext_compact_list_log2int  (send_info, &
                extra_search, mask_array, mask_shape,      &
                grid_valid_shape, src_mask, src_size, ierror)
         if (ierror /= 0) return

      else

         id_src_mask = 0
         src_mask => dummy_mask

      endif

#ifdef DEBUGX
      if ( use_source_cell ) then
         nn = 0 
         do n = 1, ncpl
            write( * ,'(a,i8)') ' Target cell ', n
            do i = 1, tgt_corners(1)
               write( * ,'(i6,3f8.3)' ) i, tgt_cell(1)%vector((i-1)*ncpl+n), &
                                           tgt_cell(2)%vector((i-1)*ncpl+n), &
                                           tgt_cell(3)%vector((i-1)*ncpl+n)
            enddo
            write( * ,'(a)') ' -------------------------- '
            do j = 1, nbr_cells(n)
               do i = 1, tgt_corners(1)
                  n1 = neighcells(nn+j,i)
                  write( * ,'(i6, 3f8.3)' ) n1, src_cell(1)%vector(n1), &
                                                src_cell(2)%vector(n1), &
                                                src_cell(3)%vector(n1)
               enddo
               if ( src_mask_available ) then
               n2 = neighcells(nn+j,1)
               write( * ,'(a,l1)' ) 'Mask set to ', src_mask(n2)
               endif
               write( * ,'(a)' ) ' ......................... '

            enddo
            nn = nn + nbr_cells(n)
            write( * ,'(a)' ) ' ========================== '
         enddo
      endif
#endif /* DEBUGX */

!
!===> Allocate vectors required for PSMILe_Trs_set_src_epio3d_dble
!
!  Create Datatype !!
!
      if ( .not. use_source_cell ) then

         do i = 1, ndim_3d
         Allocate (src_grid(i)%vector(src_size), stat = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = src_size

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'src_grid(i)%vector', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif
         end do
!
!===> Generate 3d vectors required
!
         do i = 1, ndim_2d
            call psmile_ext_compact_irreg2_dble (         &
                 send_info,                               &
                 mp%coords_pointer%coords_dble(i)%vector, &
                 mp%coords_pointer%coords_shape,          &
                 grid_valid_shape,                        &
                 src_grid(i)%vector, src_size,            &
                 ierror)

            if (ierror /= 0) return
         end do
!
         if (send_info%send_entire_valid_shape) then
!
            i = 3
            shift = grid_valid_shape (1,i) - mp%coords_pointer%coords_shape(1,i)
            len = (grid_valid_shape (2,1) - grid_valid_shape (1,1) + 1) * &
                  (grid_valid_shape (2,2) - grid_valid_shape (1,2) + 1)

            do n = 1, local_src_size / len
               src_grid(i)%vector((n-1)*len+1:n*len) = &
                    mp%coords_pointer%coords_dble(i)%vector(shift+n)
            end do

         else

            list_entries => send_info%list_entries

#ifdef PRISM_ASSERTION
            if (.not. Associated (list_entries)) then
              call psmile_assert (__FILE__, __LINE__, &
                                "list_entries are not available")
            endif
#endif

            i = 3
            shift = mp%coords_pointer%coords_shape(1,i) - 1

            src_grid(i)%vector(1:local_src_size) = &
               mp%coords_pointer%coords_dble(i)%vector( &
                             list_entries(i, 1:local_src_size) - shift)
         endif
!
!===> Add data of global search
!
         if (send_info%num2recv > 0) then
            do i = 1, ndim_3d
               ip = local_src_size
!
                  do n = 1, send_info%nrecv
                  src_grid(i)%vector(ip+1:ip+send_info%len_sent(n)) = &
                     extra_search%dble_bufs(n)%vector(           &
                                  (i-1)*send_info%len_sent(n)+1: &
                                   i   *send_info%len_sent(n))

                  ip = ip + send_info%len_sent (n)
                  end do
            end do
!
            do n = 1, send_info%nrecv
            Deallocate (extra_search%dble_bufs(n)%vector)
            end do
         endif
!
#ifdef DEBUG_TRACE
         if (ictl <= ncpl .and. use_target_grid) then
! should be searched
            ipart = 1

            if (search%grid_type == PRISM_Reglonlatvrt .or. &
                search%grid_type == PRISM_Irrlonlat_regvrt) then
               call psmile_print_irreg2_coord_dble (coords(1,ipart)%vector, &
                                                    coords(2,ipart)%vector, &
                                                    coords(3,ipart)%vector, &
                                                    shape (:, :, ipart),    &
                                                    dstijk (:, ictl), 1, 'searched')
            else if (search%grid_type == PRISM_Irrlonlatvrt) then
               call psmile_print_3d_coord_dble (coords(1,ipart)%vector, &
                                                coords(2,ipart)%vector, &
                                                coords(3,ipart)%vector, &
                                                shape (:, :, ipart),    &
                                                dstijk (:, ictl), 1, 'searched')
            endif
!
            print *, 'coords (ictl)'
            do i = 1, nb_neighbors
            if (neighbors (ictl, i) > 0) then
               print *, 'i', i, &
                  src_grid(1)%vector(neighbors (ictl, i)),&
                  src_grid(2)%vector(neighbors (ictl, i)),&
                  src_grid(3)%vector(neighbors (ictl, i))
            endif
            end do
         endif
#endif
      endif ! .not. use_source_cell

      interpolation(1:ndim_2d) = interpolation_methods(1:ndim_2d)
      interpolation(ndim_3d) = PSMILe_Undef

      call psmile_get_epio_handle (field%comp_id, grid_id, method_id, mask_id, &
                interpolation, search%msg_intersections, id_trans_out,         &
                id_trans_in, search%sender, cpl_id, epio_id, trs_rank, ierror)

      if ( epio_id /= PSMILe_undef ) then
!
!----------------------------------------------------------------------------
!  Link epio_id and id's of transient in and out if grids have already
!  been processed before.
!----------------------------------------------------------------------------
!
         call psmile_trs_set_triple_links (id_trans_out, id_trans_in, &
              epio_id, trs_rank, ierror)
         if (ierror /= 0) return

! ===> ... free memory

         if ( use_source_cell ) then
            do i = 1, ndim_3d
               Deallocate (src_cell(i)%vector)
            end do
         else
            do i = 1, ndim_3d
               Deallocate (src_grid(i)%vector)
            end do
         endif

         if (id_src_mask > 0) then
            Deallocate (src_mask)
         endif

         if ( use_target_cell ) then
            do i = 1, ndim_3d
               Deallocate (tgt_cell(i)%vector)
            end do
         else
            if (use_target_grid) then
               do i = 1, ndim_3d
                  Deallocate (tgt_grid(i)%vector)
               end do
            endif
         endif

      else
!
!----------------------------------------------------------------------------
!  Transmit complete grid information
!----------------------------------------------------------------------------
!
!===> ... Transfer of source info
!
      if ( use_source_cell ) then
         call psmile_trs_set_src_epio3d_dble (epio_id, &
              trs_rank, src_size, src_cell_size,       &
              src_cell(1)%vector,                      &
              src_cell(2)%vector,                      &
              src_cell(3)%vector,                      &
              id_src_mask, src_mask, ierror)
!
!===> ...... Free memory allocated for transfer of cells
!
            do i = 1, ndim_3d
            Deallocate (src_cell(i)%vector)
            end do
!
      else

         call psmile_trs_set_src_epio3d_dble (epio_id, &
              trs_rank, src_size, src_size,            &
              src_grid(1)%vector,                      &
              src_grid(2)%vector,                      &
              src_grid(3)%vector,                      &
              id_src_mask, src_mask, ierror)
!
!===> ...... Free memory allocated for transfer of grid points
!
            do i = 1, ndim_3d
            Deallocate (src_grid(i)%vector)
            end do
      endif

      if (ierror /= 0) return
!
!===> ... Free memory allocated for PSMILe_Trs_set_src_epio3d_dble
!
      if (id_src_mask > 0) then
         Deallocate (src_mask)
      endif
!
!----------------------------------------------------------------------------
!     Transfer data on target
!----------------------------------------------------------------------------
!
!  id_tgt_mask > 0: target mask is transferred (only for conservative)
!              = 0: target_mask is not transferred
!
      id_tgt_mask = 0

      do j = 1, n_intmethods
         if ( (interpolation_methods(j) == PSMILe_conserv2D   .or.   &
               interpolation_methods(j) == PSMILe_conserv3D ) .and.  &
               associated(search%search_mask) ) id_tgt_mask = 0
      enddo
!
      if ( id_tgt_mask == 1 ) then

         Allocate (tgt_mask(ncpl), stat = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = ncpl

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'tgt_mask', &
                 ierrp, 2, __FILE__, __LINE__ )
            return
         endif
         tgt_mask = 0
         nprev = 0
         do ipart = 1, search%npart
            if (len_cpl (ipart) > 0) then
               ibeg = nprev + 1
               iend = nprev + len_cpl (ipart)
!
               leni   = shape(2,1,ipart)-shape(1,1,ipart)
               lenij  = leni * (shape(2,2,ipart)-shape(1,2,ipart))

               do j = ibeg, iend

                  if ( dstijk(1,j) /= PSMILe_Undef .and. &
                       dstijk(2,j) /= PSMILe_Undef .and. &
                       dstijk(3,j) /= PSMILe_Undef ) then

                     n = (dstijk(3,j)-shape(1,3,ipart))*lenij + &
                         (dstijk(2,j)-shape(1,2,ipart))*leni  + &
                          dstijk(1,j)-shape(1,1,ipart)+1

                     if ( search%search_mask(ipart)%vector(n) ) tgt_mask(j) = 1

                  endif

               enddo

               nprev = nprev + len_cpl (ipart)
            endif
         end do ! ipart

      else

         id_tgt_mask = 0
         tgt_mask => dummy_mask

      endif
!
!===> ... Transfer of target info
!
      if ( use_target_cell ) then

!     remove redundant information from the target list

         nb_corners = 4
         call psmile_trs_set_tgt_epio3d_dble (                            &
              epio_id, trs_rank, ncpl, nb_corners,                        &
              tgt_cell(1)%vector, tgt_cell(2)%vector, tgt_cell(3)%vector, &
              id_tgt_mask, tgt_mask, ierror)
!
!===> ...... Free memory allocated for transfer of cells
!
            do i = 1, ndim_3d
            Deallocate (tgt_cell(i)%vector)
            end do
      else

         nb_corners = 1
         call psmile_trs_set_tgt_epio3d_dble (                            &
              epio_id, trs_rank, ncpl, nb_corners,                        &
              tgt_grid(1)%vector, tgt_grid(2)%vector, tgt_grid(3)%vector, &
              id_tgt_mask, tgt_mask, ierror)
!
!===> ...... Free memory allocated for transfer of grid points
!
         if (use_target_grid) then
            do i = 1, ndim_3d
            Deallocate (tgt_grid(i)%vector)
            end do
         endif
      endif
!
!===> ... Free memory allocated for PSMILe_Trs_set_src_epio3d_dble
!
      if ( id_tgt_mask == 1 ) Deallocate(tgt_mask)
!
!----------------------------------------------------------------------------
!     Link epio_id and id's of transient in and out
!----------------------------------------------------------------------------
!
      call psmile_trs_set_triple_links (id_trans_out, id_trans_in, &
                                        epio_id, trs_rank, ierror)
      if (ierror /= 0) return
!
!     Transfer neighbours to the transformer
!
      if ( use_source_cell ) then
         call psmile_trs_give_neighcells3d (epio_id, trs_rank,     &
              ncpl, nbr_cells, nbr_cells_tot, nb_neighbors,        &
              source_cell_index, neighcells, grid_info%grid_type,  &
              ierror)
      else
         call psmile_trs_give_neighbors3d (epio_id, trs_rank,      &
              ncpl, nb_neighbors, neighbors, ierror)
      endif
      if (ierror /= 0) return
!
!   Store Id of EPIO and of the transient to be put
!
      cpl_list(cpl_id)%epio_id  = epio_id
      cpl_list(cpl_id)%trs_rank = trs_rank

      endif ! epio_id

      send_info%epio_id  = epio_id
      send_info%trs_rank = trs_rank
!
#ifdef DEBUG
      print * , trim(ch_id), ': psmile_info_trs_loc_irreg2_dble: index, epio_id, n_list', &
               send_index, epio_id, send_info%n_list, &
               send_info%send_entire_valid_shape
#endif
!
!===> Free source locations (not needed for transfer to coupler) any more
!
1000  continue
!
      call psmile_locations_dealloc (send_info, ierror)

      if ( use_source_cell ) then
         Deallocate (source_cell_index, neighcells, nbr_cells)
      else
         Deallocate (neighbors)
      endif
!
1100  continue
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

#ifdef VERBOSE

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

9970 format (1x, a, ': psmile_info_trs_loc_irreg2_dble: eof src_size=0, ierror =', i3)

#endif /* VERBOSE */

#ifdef DEBUG
9960 format (1x, 'Taskout%In_channel(', i3, '): origin_type', i5, &
                 ', remote comp id', i4, &
                 ', global_transi_id', i4, ', remote_transi_id', i4)
#endif

      end subroutine PSMILe_info_trs_loc_irreg2_dble