psmile_find_intersect.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_Find_Intersect
!
! !INTERFACE:

      subroutine psmile_find_intersect (comp_info, global_index,     &
                                        num_intersect_per_grid,      &
                                        num_dummy_intersect_per_grid,&
                                        ninter, nmyint, nnull,       &
                                        lastag, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use psmile_grid, only : psmile_transform_extent_back, &
                              code_no_trans
!
      use PSMILe, dummy_interface => PSMILe_Find_intersect
      use psmile_common, only : enddef_msg_intersections

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

!     Index of current comp_info in all_comp_infos (:)
!
      Integer, Intent (In)               :: lastag

      Integer, Intent (InOut)            :: num_intersect_per_grid(:), 
                                            num_dummy_intersect_per_grid(:)

!     Message tag to be used
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Enddef_comp), Intent (InOut) :: comp_info

!     Info on the component (located on the actual process)
!     for which the intersections should be found.
!
      Integer, Intent (InOut)           :: ninter

!     Contains total number of intersections found for the
!     actual process.

      Integer, Intent (InOut)           :: nmyint

!     Contains total number of intersections found which are located
!     also on the actual process.

      Integer, Intent (InOut)           :: nnull

!     Contains total number of dummy (empty) intersections found
!     which were sent to a destination process
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)             :: ierror

!     Returns the error code of PSMILe_Find_intersect;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!  ... for components

      integer                         :: local_comp_id
!     local comp id of component currently processed by psmile_find_intersect
      integer                         :: global_comp_id
!     global comp id of component currently processed by psmile_find_intersect
      integer                         :: rank
!     rank of this process in comm_psmile
      integer                         :: icomp
!     loop variable for processing all collected component information
      type(enddef_comp), pointer      :: icomp_info
!     pointer to component information currently checked agains comp_info

!  ... for grids
!  ... Note: mg_gen and range are currently dimensioned with
!            Number_of_Grids_allocated. This could be decreased to the
!            range of grid id's which are used for coupling.

      integer                         :: grid_id, grid_type, extent_id
      integer                         :: inum_total_extents
!     number of extents in the currently checked component
      integer                         :: local_extents_range(2)
      integer                         :: local_grid_extents_range(2)
      integer                         :: jglob, jbeg, jend, n
      integer                         :: max_num_grids_per_comp
      integer                         :: global_dest_grid_id

      logical                         :: mg_gen (Number_of_Grids_allocated)

!  ... for methods

      integer                         :: m, method_id, n_meths
      integer                         :: poss_meths (Number_of_Methods_allocated)

!  ... for fields

      Integer                         :: ibeg, n_in_fields, n_out_fields
      Integer                         :: var_id, n_vars
      Integer                         :: poss_fields (Number_of_Fields_allocated)
      Type (ch_ptr), pointer          :: In_channel(:)

!  ... for masks

      Integer                         :: defined, last_maskid

!  ... for intersections

      Integer                         :: nbr_blocks
      Integer                         :: extent
      Integer                         :: idim(ndim_3d)
      Integer                         :: offset(ndim_3d)
      Integer                         :: addoffset(ndim_3d)
      Integer                         :: i, j, jj, jglob0
      Integer                         :: nb, ng, nl, npart, npartg
      Integer                         :: num_intersect, num_intersect_allocated
      Double Precision                :: transformed (2, ndim_3d)
      Integer                         :: range (2, ndim_3d, Number_of_Grids_allocated)
      Double Precision, Allocatable   :: dinter (:, :, :)
      Integer, Allocatable            :: inters      (:, :, :)
      Integer, Allocatable            :: inters_save (:, :, :)
      Integer, Allocatable            :: ids (:, :)
      Integer, Allocatable            :: found (:)
      Integer, Allocatable            :: buffer (:)
!
      Integer                        :: field_list (nd_field_list, 
                                                  Number_of_Fields_allocated)

      Integer                         :: field_tmp  (nd_field_list, 
                                                  Number_of_Fields_allocated)
      Integer                         :: expand

!  ... for messages containing intersections

      Integer                         :: ip
      Integer                         :: msgint  (nd_msgint)
      type (enddef_msg_intersections) :: msg_intersections

!  ... for communication

      Integer                         :: dest

!  ... for error parameters

      integer, parameter              :: nerrp = 3
      integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Find_intersect" looks for the intersections
! of the own grids with the grids of other components and
! sends the partner data on the extension of the intersection.
!
! Doxygen documentation (requires mscgen tool)
!> \file psmile_find_intersect.F90
!! General overview on tasks performed and messages sent in this routine:\n
!! \msc
!!   # entities
!!   Target,Source;
!!   ||| ;
!!   Target rbox Source [label="compute extents"] ;
!!   Target rbox Source [label="exchange extents"] ;
!!   Target rbox Target [label="find intersection"] ;
!!   Target=>Source [label = "send intersection (lasttag(100) reqidx(2))"] ;
!! \endmsc
!<
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.01.04    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_find_intersect.F90 2855 2011-01-05 17:19:54Z coquart $
! $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_find_intersect.F90 2855 2011-01-05 17:19:54Z coquart $'
!
!----------------------------------------------------------------------
!
!WARNING: Workaround
!
      double precision             :: tol = 1d-6 ! needs to be an argument

#ifdef VERBOSE
      print 9990, trim(ch_id), comp_info%comp_id, comp_info%global_comp_id

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Initialization
!
      ierror = 0
      call psmile_init_enddef_msg_inters (msg_intersections)
!
      local_comp_id  = comp_info%comp_id
      global_comp_id = comp_info%global_comp_id
      rank           = Comps(local_comp_id)%rank
!
!     local_extents_range: index range for comp_info%all_extents and
!                          comp_info%all_extent_infos in which the extents
!                          of the local process can be found
!
      local_extents_range(1) = SUM (comp_info%Number_of_grids_vector(1:rank)) + 1
      local_extents_range(2) = local_extents_range(1) + comp_info%Number_of_grids_vector(rank+1) - 1
!
!      mg_gen : Does multigrid sequence have to be generated for grid "grid_id" ?
!               mg_gen  = .false. : Multigrid sequence need not be generated
!              ? ndlev (igrid) > 0 ? Number of possible grids
!
      mg_gen (:) = .false.
!
      num_intersect_allocated = 0
!
!===> Get maximal number of extents per component for all collected components
!
      max_num_grids_per_comp = 0
      do icomp = 1, Number_of_coll_comps
         max_num_grids_per_comp = max (max_num_grids_per_comp, &
             SUM(all_comp_infos(icomp)%Number_of_grids_vector(1:all_comp_infos(icomp)%size)) )
      end do
!
!===> Allocate found vector
!
      Allocate (found (1:max_num_grids_per_comp), stat = ierror)
      if (ierror > 0) then

         call psmile_error ( PRISM_Error_Alloc, 'found', &
                             (/ierror, max_num_grids_per_comp/), &
                             2, __FILE__, __LINE__ )
         return
      endif
!
!===> For all collected components
!
      do icomp = 1, Number_of_coll_comps

         icomp_info => all_comp_infos(icomp)

#ifdef DEBUG
         print *, "testing comp: ", icomp, " of ", Number_of_coll_comps, &
                  " (global comp_id ", icomp_info%global_comp_id, ")"
#endif /* DEBUG */
         ! if the currently checked component information matches comp_info or
         ! if both components are not coupled, then we can skip this component
         if (global_comp_id == icomp_info%global_comp_id .or. &
             .not. psmile_to_be_coupled (global_comp_id, icomp_info%global_comp_id) ) cycle
!
!-----------------------------------------------------------------------
!       Loop over all parts of grids
! Tol ? um all_extents
!-----------------------------------------------------------------------
!
!  inum_total_extents = Total number of extents in component "icomp"
!
         inum_total_extents = SUM(icomp_info%Number_of_grids_vector(1:icomp_info%size))
!
         ! get the start index of the first grid of the current component
         local_grid_extents_range(1) = local_extents_range(1)
!
!        for all grids
         do while (local_grid_extents_range(1) <= local_extents_range(2))
            grid_id = comp_info%all_extent_infos(1, local_grid_extents_range(1))
            grid_type = Grids(grid_id)%grid_type
!
#ifdef PRISM_ASSERTION
            if (grid_type /= comp_info%all_extent_infos(2,local_grid_extents_range(1))) then
               print *, "grid_id, grid_type, comp_info%all_extent_infos", &
                         grid_id, grid_type, comp_info%all_extent_infos(2,local_grid_extents_range(1))
               call psmile_assert (__FILE__, __LINE__, &
                                   "Inconsistent grid type for grid grid_id")
            endif
#endif
!
!===> ... Get index range of local extents of this grid
!
            local_grid_extents_range(2) = local_grid_extents_range(1)
            do while (local_grid_extents_range(2) < local_extents_range(2))
               if (comp_info%all_extent_infos(1,local_grid_extents_range(2)+1) /= grid_id) exit
               local_grid_extents_range(2) = local_grid_extents_range(2) + 1
            enddo
!
! ... Look for an intersection of currently processed grid of actual component
!                  with a grid of icomp-th component
!
            found (1:inum_total_extents) = 0
!
!             do i = local_grid_extents_range(1), local_grid_extents_range(2)
            do i = local_grid_extents_range(1), local_grid_extents_range(2)
!cdir vector
               ! for all extents of the other component (prism_gridless
               ! grids can only be coupled among themselves)
               do j = 1, inum_total_extents
                  if (grid_types_are_compatible(grid_type,     &
                      icomp_info%all_extent_infos(2,j))  .and. &
                       comp_info%all_extents(1, 1, i) <=       &
                      icomp_info%all_extents(2, 1, j)    .and. &
                       comp_info%all_extents(2, 1, i) >=       &
                      icomp_info%all_extents(1, 1, j)    .and. &
                       comp_info%all_extents(1, 2, i) <=       &
                      icomp_info%all_extents(2, 2, j)    .and. &
                       comp_info%all_extents(2, 2, i) >=       &
                      icomp_info%all_extents(1, 2, j)    .and. &
                       comp_info%all_extents(1, 3, i) <=       &
                      icomp_info%all_extents(2, 3, j)    .and. &
                       comp_info%all_extents(2, 3, i) >=       &
                      icomp_info%all_extents(1, 3, j)) then
                     found (j) = found (j) + 1
#ifdef DEBUG
                     print *, "extents intersect"
                     print "(a,6(' ',e13.6))", " local extent :", &
                        comp_info%all_extents(:,:,i)
                     print "(a,6(' ',e13.6))", " remote extent:", &
                        icomp_info%all_extents(:,:,j)
#endif
                  endif
               end do ! j
            end do ! i
#ifdef DEBUG
            print *, "-----------------------------------------------------"
#endif
!
!===> Get total number of intersections
!
            num_intersect = sum (found (1:inum_total_extents))
            if (num_intersect == 0) then
!              ... no intersection found for grids local_grid_extents_range(1:2)
!              => set start index of next grid range
               local_grid_extents_range(1) = local_grid_extents_range(2) + 1
               cycle ! do while (local_grid_extents_range(1) <= local_extents_range(2))
            endif
!
!===> Re-allocate working vectors "dinter, inters, ids, inters_save" if necessary
!
            if (num_intersect > num_intersect_allocated) then

               if (num_intersect_allocated > 0) &
                  Deallocate (dinter, inters, ids, inters_save)
!
               Allocate (dinter      (2, ndim_3d, num_intersect+1), &
                         inters      (2, ndim_3d, num_intersect),   &
                         inters_save (2, ndim_3d, num_intersect),   &
                         ids (num_intersect, 2), stat = ierror)
               if (ierror /= 0) then
                  call psmile_error (                                                &
                     PRISM_Error_Alloc, 'dinter,inters,ids, inters_save',            &
                     (/ierror, 3*(2*ndim_3d*(num_intersect+1))+num_intersect*2/), 2, &
                     __FILE__, __LINE__ )
                  return
               endif
               num_intersect_allocated = num_intersect
            endif
!
!===> ... Get list of methods id's which may be used for coupling.
!
!  n_in_fields  = Number of input  fields for this grid;
!                 i.e. number of fields which receive data from
!                      other applications.
!  n_out_fields = Number of output fields for this grid;
!                 i.e. number of fields for which data of other
!                      applications must be searched and interpolated.
!
            n_in_fields = 0
            n_out_fields = 0
            n_meths = 0
!
!cdir vector
            if (Grids(grid_id)%comp_id == local_comp_id) then
               do i = 1, Number_of_Methods_allocated
                  if (Methods(i)%grid_id == grid_id               .and. &
                     Methods(i)%status  == PSMILe_status_defined .and. &
                     Methods(i)%used_for_coupling ) then
                     n_meths = n_meths + 1
                     poss_meths (n_meths) = i
                  endif
               end do
            endif
!
            if (n_meths == 0) then
!                 Grid is not used for coupling
!                 test before !!
!                 Not so easy, as the Grid can be used as Taskout grid only
#ifdef DEBUG
               print *, "n_meths ", n_meths
               write (*, 9960) trim(ch_id), local_comp_id, grid_id
#endif /* DEBUG */
            endif
!
!===> ... Get list of input fields of grid "grid_id"
!         which may be used for coupling
!         (*) Look for a input field of component "local_comp_id" and grid "grid_id"
!             which is used for coupling
!         (*) Control whether this input field has to get data from component
!             "icomp_info%local_comp_id"
!             TODO: ?? Could we also check the Grid ID ?
!                   ?? This is not contained in structure ch_ptr ! ?
!
            if (n_meths > 0) then
!cdir vector
               do i = 1, Number_of_Fields_allocated
                  if (Fields(i)%status  == PSMILe_status_defined) then
                     if ( Fields(i)%comp_id == local_comp_id                    .and. &
                        Methods(Fields(i)%method_id)%grid_id == grid_id .and. &
                        Fields(i)%used_for_coupling ) then

                        IF ( ASSOCIATED(Fields(i)%Taskin%In_channel) ) THEN
                            DO j = 1, SIZE(Fields(i)%Taskin%In_channel)
                              IF (Fields(i)%Taskin%In_channel(j)%remote_comp_id == &
                                 all_comp_infos(icomp)%global_comp_id ) THEN
                                  n_in_fields = n_in_fields + 1
                                  poss_fields (n_in_fields) = i
                              ENDIF
                            ENDDO
                        ENDIF
!===>   Get number of output fields
                        if ( Associated(Fields(i)%Taskout) ) then
                            DO j = 1, SIZE(Fields(i)%Taskout)
                              IF (Fields(i)%Taskout(j)%remote_comp_id == &
                                 all_comp_infos(icomp)%global_comp_id ) THEN
                                  n_out_fields = n_out_fields + 1
                              ENDIF
                            ENDDO
                        endif

                     endif
                  endif
               enddo
!
#ifdef DEBUG
               if (n_in_fields + n_out_fields == 0) then
!                 Grid is not used for coupling
!                 test before !!
                  print *, "n_in_fields, n_out_fields ", n_in_fields, n_out_fields
                  write (*, 9960) trim(ch_id), local_comp_id, grid_id
               endif
#endif /* DEBUG */
            endif
!
!===>
!
            jbeg = 1

            do while ( jbeg <= inum_total_extents )
!
               do jglob = jbeg, inum_total_extents
                  if (found (jglob) > 0) exit
               end do
!
               if (jglob > inum_total_extents) exit
!
!===> ... non empty intersection between local grids "local_grid_extents_range(1:2)" and
!                                        global grid "jglob" found
!
!===> ... Get destination process in communicator comm_psmile
!         !!! We need something better !!!
!
               n = 0
!cdir vector
               do i = 1, icomp_info%size
                  if (jglob <= n+icomp_info%Number_of_grids_vector(i)) exit
                  n = n + icomp_info%Number_of_grids_vector(i)
               end do

               if (i > icomp_info%size) then
                  ierror = PRISM_Error_Internal
                  ierrp (1) = icomp
                  ierrp (2) = icomp_info%size
                  call psmile_error (ierror, "Cannot find grid", ierrp, 2, &
                                     __FILE__, __LINE__)
                  return
               endif
!
               dest = icomp_info%psmile_ranks(i)
               jend = n + icomp_info%Number_of_grids_vector(i)
!
!===> ... Get number of (logical) partitions of global grid "jglob"
!         (cf. psmile_enddef_comp.F90)
!         all_extent_infos (1, *) = Local grid id (extent id)
!
               extent_id = icomp_info%all_extent_infos(1,jglob)
!
               do jbeg = jglob+1, jend
                  if (icomp_info%all_extent_infos(1,jbeg) /= &
                      extent_id) exit
               end do
!
               if (n_in_fields + n_out_fields == 0) then
                  npart = 0
                  go to 1000
               endif
!
               npartg = jbeg - jglob
               jglob0 = jglob - 1
!
               global_dest_grid_id = icomp_info%all_extent_infos(4,jglob)
!
!===> ... Determine intersections
! minval > 0 leads to problems for 2d intersections
! therefore minval >= 0
!
               n = 0
               do ng = 1, npartg
                  if (found (jglob0+ng) > 0) then
                     do i = local_grid_extents_range(1), local_grid_extents_range(2)
                        dinter (1,:, n+1) = max ( &
                               comp_info%all_extents(1,:,i), &
                              icomp_info%all_extents(1,:,jglob0+ng))
                        dinter (2,:, n+1) = min ( &
                               comp_info%all_extents(2,:,i), &
                              icomp_info%all_extents(2,:,jglob0+ng))

                        if (minval (dinter (2,:, n+1) - dinter (1,:, n+1)) >= 0.0d0) then
!
                           n = n + 1
                           ids (n, 1) = i
                           ids (n, 2) = jglob0 + ng
#ifdef DEBUG
                           print "(a,6(' ',e13.6))", "local extent :", &
                              comp_info%all_extents(:,:,i)
                           print "(a,6(' ',e13.6))", "remote extent:", &
                              icomp_info%all_extents(:,:,jglob0+ng)
                           print "(a,6(e13.6,' '))", "intersection: ", dinter (:,:, n)
#endif
                        endif
                     end do !  i = local_grid_extents_range(1), local_grid_extents_range(2)
                  endif
               end do ! npartg
!
#ifdef PRISM_ASSERTION
               if (n > num_intersect) then
                  print *, 'n, num_intersect', n, num_intersect
                  call psmile_assert (__FILE__, __LINE__, 'n > num_intersect !')
               endif
#endif
               npart = n
!
!===> ...  Transfer intersection of grid "grid_id" (igrid) back if necessary.
!
               do n = 1, npart
                  if (comp_info%all_extent_infos(3, ids(n,1)) /= code_no_trans) then
                     call psmile_transform_extent_back ( &
                           (/comp_info%all_extent_infos(3,ids(n,1))/), &
                           dinter(:, :, n), transformed, 1, ierror)
                     if ( ierror /= 0 ) return
                     dinter (:, :, n)  = transformed
#ifdef DEBUG
                     print "(a,6(e13.6,' '),a)", "intersection: ", dinter (:,:, n), " (transformed)"
#endif
                     endif
               end do
!
!===> ... Determine sub range of intersection "dinter" using grid coordinates
!
!  ?!? What about 8 subregions, and if we only send those subregions,
!  which have an intersection ?!?!
!
               if (grid_type == PRISM_Gridless) then
!
!===> ........... Gridless grid
!                 Note: Initially, the send and receive areas should
!                 be locally at this place. Since masks can be
!                 assigned to gridless grids, the data on gridless
!                 grids is sent as usual to the counterpart.
!                 TODO: Send data only if masks are used on any process.
!
                  inters (:, :, 1:npart) = nint (dinter (:, :, 1:npart))

                  inters_save = inters
!
!                 Transform intersection into local indices (move into a new subroutine)
!                 ToDo: move into a new subroutine:
!
!                     functionality required by
!                     - psmile_find_intersect
!                     - psmile_search_donor_gridless
!                     - psmile_get_locations_3d
!
                  if (Associated (Grids(grid_id)%partition)) then

                     nbr_blocks = size(Grids(grid_id)%partition(:,1))

                     if ( nbr_blocks > 1 ) then

                        do n = 1, npart

                           ! find correct block

                           do nb = 1, nbr_blocks

                              if ( inters(2,1,n) >= Grids(grid_id)%partition(nb,1) + 1 .and. &
                                   inters(1,1,n) <= Grids(grid_id)%partition(nb,1) +         &
                                                    Grids(grid_id)%extent   (nb,1)     .and. &
                                   inters(2,2,n) >= Grids(grid_id)%partition(nb,2) + 1 .and. &
                                   inters(1,2,n) <= Grids(grid_id)%partition(nb,2) +         &
                                                    Grids(grid_id)%extent   (nb,2)     .and. &
                                   inters(2,3,n) >= Grids(grid_id)%partition(nb,3) + 1 .and. &
                                   inters(1,3,n) <= Grids(grid_id)%partition(nb,3) +         &
                                                    Grids(grid_id)%extent   (nb,3) ) exit

                           enddo

                           ! transform into local indices

                           if ( nb <= nbr_blocks ) then
#ifdef DEBUG
                              print *, ' Transform block ', nb, Grids(grid_id)%partition(nb,1)+1 &
                                                   , ' - ',     Grids(grid_id)%partition(nb,1)+  &
                                                                Grids(grid_id)%extent   (nb,1)
                              print *, ' global inter    ', nb, inters(:, :, n)
#endif
                              offset    = 0
                              addoffset = 0

                              do i = 1, ndim_3d
                                 idim(i) = Grids(grid_id)%grid_shape(2,i) - &
                                           Grids(grid_id)%grid_shape(1,i)+1
                                 do nl = 1, nb-1
                                    addoffset(i) = Grids(grid_id)%extent(nl,i) + offset(i)
                                    if ( addoffset(i) < idim(i) ) &
                                         offset(i) = addoffset(i)
                                 enddo
                              enddo

                              do i = 1, ndim_3d
                                 extent = inters(2,i,n) - inters(1,i,n)
                                 inters(1, i, n) = inters(1, i, n) + offset(i) &
                                                 - Grids(grid_id)%partition(nb,i) &
                                                 + Grids(grid_id)%grid_shape(1,i) - 1
                                 inters(2, i, n) = inters(1, i, n) + extent
                              end do
#ifdef DEBUG
                              print *, ' part offset ', nb, offset(:)
                              print *, ' local inter ', nb, inters(:, :, n)
#endif
                           else

                              ierror = PRISM_Error_Internal
                              ierrp (1) = nb
                              ierrp (2) = nbr_blocks
                              call psmile_error (ierror, "No block found", ierrp, 2, &
                                   __FILE__, __LINE__)
                              return

                           endif

                        enddo

                     else

                        ! case for single block

                        do n = 1, npart
                           do i = 1, ndim_3d
                              inters (1:2, i, n) = &
                              inters (1:2, i, n) - Grids(grid_id)%partition (1,i) &
                                                 + Grids(grid_id)%grid_shape(1,i) - 1
                           end do
                        end do

                     endif

                  endif
!
                  if (npart > 1) then
                     call psmile_remove_intersect_int (inters, ids(1:npart,1), &
                                                       ids(1:npart,2), npart, ierror)
                     if (ierror > 0) return
                  endif

#ifdef PRISM_ASSERTION
                  if (icomp_info%all_extent_infos(2,jglob) /= &
                      PRISM_Gridless) then
                     call psmile_assert (__FILE__, __LINE__, &
                      'Coupling of gridless grid with real grid is impossible !')
                  endif
#endif
               else
!
!===> ........... Real grid
!                 Remove common intersections
!
                  if (npart > 1) then
                     call psmile_remove_intersect (dinter, ids(1:npart,1), &
                               ids(1:npart,2), npart,                      &
                               comp_info%all_extent_infos,                 &
                               icomp_info%all_extent_infos,                &
                               ierror)
                     if (ierror > 0) return
                  endif
!
                  do n = 1, npart
#ifdef DEBUG
                       PRINT *, ' In psmile_find_intersect intersection :',n,'over :',npart
#endif
                     call psmile_sel_grid_range (grid_id, dinter (:, :, n), &
                                                 inters(:, :, n), ierror)
                     if ( ierror /= 0 ) return
                  end do
!
                  n = 1
                  do while (n <= npart)
                     if (minval (inters(2, :, n) - inters (1, :, n)) < 0) then
                        if (n /= npart) then
                           inters (:, :, n) = inters (:, :, npart)
                           ids (n, 1) = ids (npart, 1)
                           ids (n, 2) = ids (npart, 2)
                        endif

                        npart = npart - 1
                     else
                        n = n + 1
                     endif
                  end do ! while (n <= nparts)
#ifdef HUHU
!
!===> ... Remove common intersections in "inters"
!         !!! Note: This was done for the reduced gaussian grid.
!         !!! TODO
!         !!! The better (correct) solution would be to generate
!         !!! a list of points to be searched, which should be
!         !!! created here in PSMILe_sel_grid_range and which
!         !!! should be managed locally if possible.
!         !!! Transformation code for list of point and removed
!             intersections !!!
!
                  if (npart > 1) then
                     call psmile_remove_intersect_int (inters, ids(1:npart,1), &
                                 ids(1:npart,2), npart, ierror)
                     if (ierror > 0) return
                  endif
#endif
               endif
!
!===> ... Get list of varid's for each method in "list_of_methods"
!         (*) Look for a input field of component "local_comp_id", grid "grid_id" and
!                                       method "method_id",
!             which is used for coupling
!         (*) Control whether this input field has to get data from component
!             "icomp_info%comp_id"
!             TODO: ?? Could we also check the Grid ID ?
!                   ?? This is not contained in structure ch_ptr ! ?
!
               if (n_in_fields == 0) go to 1000
!
               n_vars = 0
!
               do m = 1, n_meths
                  ibeg = 1
                  method_id = poss_meths (m)
!
search_vars:      do while (ibeg <= n_in_fields)
!
!                    Look for a matching field
!
!cdir vector
                     do i = ibeg, n_in_fields
                        if (Fields(poss_fields(i))%method_id == method_id) exit
                     end do
!
                     ibeg = i + 1
                     if (i > n_in_fields) exit search_vars
!
!                    Is this field input field of the component
!                    "icomp_info%global_comp_id" ?
!
#ifdef DEBUG
                     print *, 'search Taskin for dest component', icomp_info%global_comp_id
#endif /* DEBUG */
                     In_channel => Fields(poss_fields(i))%Taskin%In_channel
!
                     do j = 1, size (In_channel)
#ifdef DEBUG
                        print 9970, j, &
                                    In_channel(j)%origin_type, &
                                    In_channel(j)%remote_comp_id, &
                                    In_channel(j)%global_transi_id, &
                                    In_channel(j)%remote_transi_id, &
                                    global_dest_grid_id
#endif /* DEBUG */
                        if (In_channel(j)%origin_type    == PSMILe_comp .and. &
                            In_channel(j)%remote_comp_id ==                   &
                                                   icomp_info%global_comp_id  &
                            .and.                                             &
                            field2grid(In_channel(j)%remote_transi_id) ==     &
                                        global_dest_grid_id) then
                           n_vars = n_vars + 1
                           field_tmp (1, n_vars) = method_id
                           field_tmp (2, n_vars) = poss_fields(i)
                           field_tmp (4, n_vars) = In_channel(j)%global_transi_id
                           field_tmp (5, n_vars) = In_channel(j)%remote_transi_id
                           field_tmp (6, n_vars) = 0
                           do jj = 1, 3
                             if ( In_channel(j)%interp%interp_meth(jj) == PSMILe_conserv2D .or. &
                                  In_channel(j)%interp%interp_meth(jj) == PSMILe_conserv3D ) then
                                  field_tmp (6, n_vars) = In_channel(j)%interp%interp_meth(jj)
                             endif
                           enddo
                        endif
                     end do ! j
!
                  end do search_vars ! while
               end do ! m
!rr TODO:
!rr silly sort of field_list to have conservative remapping first in the list
!rr This list needs to be revised. It should be arranged by method_ids, where
!rr for each method_id the conservative remapping is at the beginning.
!rr
               j = 0
               do n = 1, n_vars
                  if ( field_tmp (6, n) == PSMILe_conserv2D .or. &
                        field_tmp (6, n) == PSMILe_conserv3D ) then
                     j = j + 1
                     field_list(:,j) = field_tmp(:,n)
                  endif
               enddo
               do n = 1, n_vars
                  if ( field_tmp (6, n) /= PSMILe_conserv2D .and. &
                        field_tmp (6, n) /= PSMILe_conserv3D ) then
                     j = j + 1
                     field_list(:,j) = field_tmp(:,n)
                  endif
               enddo
!
               if (n_vars == 0) then
!                 Grid is not used for coupling with component icomp_info%comp_id
!                 test before !!
#ifdef DEBUG
                  print *, "n_vars", n_vars
                  write (*, 9960) trim(ch_id), local_comp_id, grid_id
#endif /* DEBUG */
!
                  npart = 0
               endif
!
!===> ... Is it necessary to generate multi grids for search of donor cells ?
!
1000           continue
!

               if (n_out_fields > 0 .and. npart > 0 .and. &
                  grid_type /= PRISM_Gridless) then
                  if (.not. mg_gen (grid_id)) then
                     mg_gen (grid_id) = .true.
                     range (:, :, grid_id) = inters (:, :, 1)
                  endif
!
!              ... Get bounding box for all intersections
!                  TODO: Something better (for example several parts)
!
                  do n = 1, npart
                     range (1, 1:ndim_3d, grid_id) = min ( &
                     range (1, 1:ndim_3d, grid_id), inters (1, 1:ndim_3d, n))
                     range (2, 1:ndim_3d, grid_id) = max ( &
                     range (2, 1:ndim_3d, grid_id), inters (2, 1:ndim_3d, n))
                  end do
               endif
!
!===> ... Dummy intersection found ?
!         Increment number of intersections

#ifdef DEBUG
               PRINT*, 'Number of input fields :',n_in_fields
               PRINT*, 'Number of output fields :',n_out_fields
               call psmile_flushstd
#endif
!
               if (npart == 0 .or. n_in_fields == 0) then
                  if (dest == psmile_rank) continue

                  nnull = nnull + 1
                  num_dummy_intersect_per_grid(grid_id) = &
                        num_dummy_intersect_per_grid(grid_id) + 1
                  n_vars = 0
                  npart = 0
               endif
!
               ninter = ninter + 1
               num_intersect_per_grid(grid_id) = &
                     num_intersect_per_grid(grid_id) + 1
#define CONS_REMAP_DEADLOCK_WORKAROUND
#ifdef CONS_REMAP_DEADLOCK_WORKAROUND
               paction%intersect_ranks(ninter) = dest
#endif
!
               if (dest == psmile_rank) nmyint = nmyint + 1
!
!===> ... Send message to the destination process
!
!===> ... Relative tag/counter
!
               msg_intersections%relative_msg_tag = ninter - nnull
!
!===> ... Info on target component and grid.
!         This is the component which has put the field data.
!
               msg_intersections%src_comp_id = icomp_info%comp_id
!
!===> ... Info on sending component and grid.
!         This is the component which get the field data.
!
               msg_intersections%tgt_comp_id = local_comp_id
               msg_intersections%tgt_grid_id = grid_id

               msg_intersections%all_comp_infos_comp_idx = global_index
!
               msg_intersections%num_parts = npart
!
               if (n_vars > 0) then
!
!         TODO: In der field_list fehlt noch der datentype der Methode !!!
                  do i = 1, n_vars
                     field_list (3, i) = Fields(field_list(2,i))%mask_id
                  end do
!
!              Control mask values in intersections
!                 defined == 2: All mask values are true
!                 defined == 0: All mask values are false
!
!              TODO: improve handling of already treated Masks (see last_maskid)
!                    (already_controlled (1:nvars) or something like that)
!              !!! For now this is only done for PRISM_Gridless, because it is
!                  not clear to me how much CPU time this additional control will cost. !!!
!rr            !rr For conservative remapping we have to send the mask only once.
!rr                if we do it twice, once when sending the corners and a second time
!rr                when sending to coords this will corrupt the memory. Therefore I
!rr                activated the check for all grids.
!rr
                  do n = 1, n_vars
                     if ( field_list(6,n) == PSMILe_conserv2D .or. &
                        field_list(6,n) == PSMILe_conserv3D ) then
                        last_maskid = field_list(3,n)
                        do i = n+1, n_vars
                           if ( field_list(3,i) == last_maskid ) then
                              field_list(3,i) = PRISM_UNDEFINED
                           endif
                        enddo
                     endif
                  enddo

                  if (grid_type == PRISM_Gridless) then
                     last_maskid = PRISM_UNDEFINED
!
                     do i = 1, n_vars
                     if (field_list (3, i) /= PRISM_UNDEFINED) then
                        if (field_list (3, i) == last_maskid) cycle

                        call psmile_is_mask_defined (                &
                           Masks(field_list (3, i))%mask_array, &
                           Masks(field_list (3, i))%mask_shape, &
                           inters, npart, defined, ierror)
                        if (ierror > 0) return

!                 ?? Skip in case defined == 0 ??
!
                        if (defined == 2) then
                              do j = i+1, n_vars
                              if (field_list (3,j) == field_list (3,i)) &
                                 field_list (3,j) = PRISM_UNDEFINED
                              end do
!
                           field_list (3,i) = PRISM_UNDEFINED
                        else
                           last_maskid = field_list (3,i)
                        endif

                     endif
                     end do
                  endif
!
!===> ...... Store data on the first field id
!
                  method_id = field_list (1, 1)
                  var_id    = field_list (2, 1)
!
                  msg_intersections%first_tgt_method_id = method_id
                  msg_intersections%first_tgt_var_id = var_id
!
                  msg_intersections%tgt_mask_id = field_list (3, 1)
!
!===> ...... Transfer transient id for Epio
!
                  msg_intersections%transient_in_id        = field_list (4, 1)
                  msg_intersections%transient_out_id       = field_list (5, 1)
                  msg_intersections%num_vars               = n_vars
                  msg_intersections%requires_conserv_remap = field_list (6, 1)
!
!===> ..... Get method values
!
                  msg_intersections%method_type = Methods (method_id)%method_type

                  select case (Methods (method_id)%method_type)
                  case (PSMILe_PointMethod)
                     msg_intersections%method_datatype = &
                           Methods (method_id)%coords_pointer%coords_datatype

                  case (PSMILe_VectorPointMethod)
#ifdef TODO
                     msg_intersections%method_datatype = &
                           Methods (method_id)%vector_pointer%vector_datatype
#else
                     ierrp (1) = Methods (method_id)%method_type
                     ierror = PRISM_Error_Internal

                     call psmile_error ( ierror, 'vector method type is currently not supported', &
                                      ierrp, 1, __FILE__, __LINE__ )
#endif

                  case (PSMILe_SubgridMethod)
                     msg_intersections%method_datatype = &
                           Methods (method_id)%subgrid_pointer%subgrid_datatype

                  case default
                     ierrp (1) = Methods (method_id)%method_type
                     ierror = PRISM_Error_Internal

                     call psmile_error ( ierror, 'unsupported method type', &
                                       ierrp, 1, __FILE__, __LINE__ )
                  end select

               endif ! (n_vars > 0)
!
!===> ... Data on intersection
!
               ip = ip_msgint_inter + npart*(2*ndim_3d)
!
! The fields to be handled are ordered such that those requireing conservative remapping
! are listed first. For conservative remapping the intersections have to be increased by
! 1. When the non-conservative remapping fields are processed we have to switch back to
! the original intersections.
!
               expand = PSMILe_Undef

               if (npart > 0) then

                  msg_intersections%first_src_all_extents_grid_id = ids (1, 2) ! jglob0 + ng
                  msg_intersections%first_tgt_all_extents_grid_id = ids (1, 1) ! local_grid_extents_range(1) + x

                  if (grid_type == PRISM_Gridless) then
                     allocate (msg_intersections%intersections(2*npart))
                  else
                     allocate (msg_intersections%intersections(npart))
                  endif
!
!              increase inters by one for exchange of corner info. For Gaussreduced grids
!              we currently send complete (redundant corner information. The arrays sizes
!              for Gaussreduced corners are adjusted in the repective routines.
!
                  if ( grid_type /=  PRISM_Gaussreduced_regvrt ) then

                     if ( msg_intersections%requires_conserv_remap == PSMILe_conserv2D ) then

                        do n = 1, npart
                           inters (2, 1:ndim_2d, n) = inters (2, 1:ndim_2d, n) + 1
                        enddo

                        expand = PSMILe_conserv2D

                     else  if ( msg_intersections%requires_conserv_remap == PSMILe_conserv3D ) then

                        do n = 1, npart
                           inters (2, 1:ndim_3d, n) = inters (2, 1:ndim_3d, n) + 1
                        enddo

                        expand = PSMILe_conserv3D

                     endif

                     ! MoHa: is this ever executed?
                     if ( expand == PSMILe_conserv2D .and. &
                          msg_intersections%requires_conserv_remap /= PSMILe_conserv2D ) then

                        do n = 1, npart
                           inters (2, 1:ndim_2d, n) = inters (2, 1:ndim_2d, n) - 1
                        enddo

                     else  if ( expand == PSMILe_conserv3D .and. &
                                msg_intersections%requires_conserv_remap /= PSMILe_conserv3D ) then

                        do n = 1, npart
                           inters (2, 1:ndim_3d, n) = inters (2, 1:ndim_3d, n) - 1
                        enddo

                     endif

                  endif ! ( grid_type /=  PRISM_Gaussreduced_regvrt )

                  do n = 1, npart
                     msg_intersections%intersections(n)%intersection = inters(:, :, n)
                     msg_intersections%intersections(n)%tgt_all_extents_grid_id = ids (n, 1)
                     msg_intersections%intersections(n)%src_all_extents_grid_id = ids (n, 2)
                  enddo
!
                  ip = ip + 2 * npart

                  if (grid_type == PRISM_Gridless) then
!
!===>             Transfer offset
!
                     if (Associated(Grids(grid_id)%partition)) then

                        do n = 1, npart

                           msg_intersections%intersections(n+npart)%intersection = inters_save(:, :, n)
                           msg_intersections%intersections(n+npart)%tgt_all_extents_grid_id = psmile_undef
                           msg_intersections%intersections(n+npart)%src_all_extents_grid_id = psmile_undef
                        enddo
                     endif

                     ip = ip + 1 + npart*2*ndim_3d

!                      if (npart > maxpart ) then
! #ifdef DEBUG
!                         print *, ' npart is ', npart , '!'
!                         print *, ' Preparing send of an additional message! '
! #endif
!                         allocate (buffer(2*ndim_3d*(npart-maxpart)), stat = ierror)
!                         if (ierror > 0) then
!                            ierrp (1) = ierror
!                            ierrp (2) = 2*ndim_3d*(npart-maxpart)
!
!                            ierror = PRISM_Error_Alloc
!
!                            call psmile_error ( ierror, 'found', &
!                               ierrp, 2, __FILE__, __LINE__ )
!                            return
!                         endif
!
!                         buffer(:) = pack (inters_save(:, :, npart-maxpart:npart), .true.)
!
!                      endif

                  endif
               else
                  msg_intersections%first_tgt_all_extents_grid_id = psmile_undef
                  msg_intersections%first_src_all_extents_grid_id = psmile_undef
                  nullify (msg_intersections%intersections)
               end if ! (npart > 0)
!
#ifdef DEBUG
               print *, "psmile_find_intersect: psmile_rank->dest, tag:", &
                        psmile_rank, dest, lastag
               PRINT *, "psmile_find_intersect: datatype, npart:", &
                        msg_intersections%method_datatype, npart
               PRINT*, 'psmile_find_intersect : ninter, nmyint , nnull:',ninter, nmyint , nnull
               call psmile_flushstd
               do n = 1, npart
                  print *, "psmile_find_intersect: ipart, inter", n, inters (:, :, n)
                  print *, "psmile_find_intersect: ipart, ids", n, ids(n, 1), ids (n,2)
                  call psmile_flushstd
               end do
#endif /* DEBUG */
!
!===> ... Send data to the process which contains intersection
!         The corresponding receive is found in psmile_enddef_action
!
#ifdef DEBUG
               print *, ' Sending tag ', lastag, ' to destination ', dest
#endif /* DEBUG */

               ! copy collected data into a sendable buffer
               call psmile_pack_msg_intersections(msg_intersections, msgint)

               call psmile_bsend (msgint, ip, MPI_INTEGER, dest, lastag, &
                                  comm_psmile, ierror)

               if (ierror /= MPI_SUCCESS) then
                  ierrp (1) = ierror
                  ierrp (2) = dest
                  ierrp (3) = lastag
                  ierror = PRISM_Error_Send

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

               if ( associated (msg_intersections%intersections)) &
                  deallocate (msg_intersections%intersections)

               if (grid_type == PRISM_Gridless .and. npart > maxpart ) then

#ifdef PRISM_ASSERTION
                  call psmile_assert (__FILE__, __LINE__, &
                                   "This part of the code does not work. I am sorry.")
#endif

!                   call psmile_bsend (buffer, 2*ndim_3d*(npart-maxpart), &
!                      MPI_INTEGER, dest, ind_msgint_tag, comm_psmile, ierror)
!                   if (ierror /= MPI_SUCCESS) then
!                      ierrp (1) = ierror
!                      ierrp (2) = dest
!                      ierrp (3) = lastag
!                      ierror = PRISM_Error_Send
!
!                      call psmile_error (ierror, 'MPI_Send', &
!                         ierrp, 3, __FILE__, __LINE__ )
!                      return
!                   endif

               endif
!
!===> ... Send data to the process containing data on further fields
!         TODO: In field_list the datatype of Method is still missing !!!
!
               if (n_vars > 1) then
#ifdef DEBUG
                  print *, ' Sending tag ', vartag, ' to destination ', dest
#endif /* DEBUG */
                  call psmile_bsend (field_list(1,2), (n_vars-1)*nd_field_list, &
                                    MPI_INTEGER, dest, vartag, comm_psmile,  &
                                    ierror)
                  if (ierror /= MPI_SUCCESS) then
                     ierrp (1) = ierror
                     ierrp (2) = dest
                     ierrp (3) = vartag
                     ierror = PRISM_Error_Send

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

            end do ! while ( jbeg <= inum_total_extents )

            ! set index of next grid range
            local_grid_extents_range(1) = local_grid_extents_range(2) + 1
         enddo ! while (local_grid_extents_range(1) < local_extents_range(2))

      end do ! all components
!
!===> Prepare search for grid with id "grid_id"
!     ? Delay until first receive ?
!
      do grid_id = 1, Number_of_Grids_allocated
         if (mg_gen (grid_id)) then
!
            call psmile_mg_setup (grid_id, range(1, 1, grid_id), &
                                 tol, ierror)
            if (ierror /= 0) return
!
         endif
      end do
!
!=======================================================================
!     Free
!=======================================================================
!
      Deallocate (found)
      if (num_intersect_allocated > 0) then
         Deallocate (dinter, inters, ids, inters_save)
      endif
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror
      call psmile_flushstd
#endif /* VERBOSE */

!
!  Formats:
!

#ifdef VERBOSE

9990 format (1x, a, ': psmile_find_intersect: local comp_id', i3, ' global comp_id', i3)
9980 format (1x, a, ': psmile_find_intersect: eof ierror =', i3)

#endif /* VERBOSE */

#ifdef DEBUG
9970 format (1x, 'Taskin%In_channel(', i3, '): origin_type', i5, &
                 ', remote comp id', i4, &
                 ', global_transi_id', i4, ', remote_transi_id', i8, &
            /1x, 'global_dest_grid_id', i4)
9960 format (1x, a, ': psmile_find_intersect: #### Grid is not used for coupling !', &
            /1x, 'Introduce test before !!! comp_id', i4, '; grid id', i4)
#endif /* DEBUG */

      contains

         ! this function needs to be moved to a central place!
         function grid_types_are_compatible(type1, type2)

            use prism, only: prism_gridless

            integer, intent(in) :: type1, type2
            logical grid_types_are_compatible

            ! gridless grids can only be coupled among themselves
            grid_types_are_compatible = &
               type1 == prism_gridless .eqv. type2 == prism_gridless

         end function grid_types_are_compatible

      end subroutine PSMILe_Find_intersect