psmile_add_points_found.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_Add_points_found
!
! !INTERFACE:

      subroutine psmile_add_points_found (grid_id, search,   &
                           extra_search, indices_returned,   &
                           found, n_send, len_nsend,         &
                           index_found, n_found,             &
                           neighbors_3d, nloc, num_neigh,    &
                           grid_valid_shape, use_how, ierror)
!
! !USES:
!
      use PRISM
!
      use PSMILe, dummy_interface => PSMILe_Add_points_found
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer,           Intent (In)     :: grid_id
!
!     Handle to the pointer to the receive list
!
      Type (Enddef_search),     Intent (InOut) :: search

!     Info's on coordinates to be searched
!
      Integer,           Intent (In)     :: n_send
!
!     Number of points for which additional interpolation data was found.
!     Dimension of "found" and "indices_returned".

      Integer,            Intent (In)    :: indices_returned (n_send)
!
!     Indices of the extra points in vectors "required" and
!     "indices_req" for which interpolation data was found.
!     This is a local index of the partition "ipart".
!     (cf. routine PSMILe_Store_faces_3d_reg_dble)

      Integer,            Intent (In)    :: n_found
!
!     Total number of points found which are used for interpolation.
!     Dimension of "index_found".
!
      Integer,            Intent (In)    :: index_found (n_found)
!
!     Index of point in compact list "liste"
!
      Integer,            Intent (In)    :: len_nsend (search%npart)
!
!     Number of points for which additional interpolation data was required
!     per partition.
!
      Integer,            Intent (In)    :: nloc
!
!     Total  Number of locations to be transferred
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.

      Integer,            Intent (In)    :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.

      Integer,            Intent (In)    :: grid_valid_shape (2, ndim_3d)

!     Specifies the valid block shape for the "ndim_3d"-dimensional block

      Integer,            Intent (In)    :: use_how

!     Information about how to apply the mask

! !INPUT/OUTPUT PARAMETERS:
!
      Type (Extra_search_info), Intent (InOut) :: extra_search
!
!     Number of locations where
!       (*) required mask values were not true

      Integer,            Intent (InOut) :: found    (n_send)
!
!     Input: Code of the interpolation bases found for corresponding point.
!            (cf. "required")
!     Output: Code cleaned by points masked out.

      Integer,            Intent (InOut) :: neighbors_3d (ndim_3d, nloc, 
                                                          num_neigh)
!
! !OUTPUT PARAMETERS:
!
      Integer,            Intent (Out)   :: ierror

!     Returns the error code of PSMILe_Add_points_found;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS:
!
! masked_out   = Hash value for a point which is masked out
!
      Integer, Parameter           :: masked_out = 0
      Integer, Parameter           :: undef      = -12
!
! !LOCAL VARIABLES
!
!     ... loop variables
!
      Integer                         :: i, k, n
!
!     ... Bit value for found points
!
      Integer                         :: iband (n_send)
      Integer                         :: code
!
!     ... 
!
      Integer, Allocatable            :: corner2ind (:, :)
      Integer                         :: n_corners
      Integer                         :: outside
      Integer, Pointer                :: recv_list(:)
!
!     ... for partitions
!
      Integer                         :: ipart, ip
      Integer, Pointer                :: indices_req (:)
      Integer, Pointer                :: required (:)
      Integer, Pointer                :: len_req (:)
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Add_points_found" add information on
! points found by global search to the already existing information.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 02.02.06    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_add_points_found.F90 2935 2011-02-03 07:52:41Z redler $
!  $Author: redler $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_add_points_found.F90 2935 2011-02-03 07:52:41Z redler $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), n_send, n_found

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      if (associated(Grids(grid_id)%recv_list)) then
          !
          ! in particular for Gauss-reduced grids
          !
          recv_list => Grids(grid_id)%recv_list
          outside   = grid_valid_shape(2,1) + sum(recv_list) + 11
      else
          outside   = grid_valid_shape(2,1) + 11
      endif

      n_corners = num_neigh
!
      len_req     => extra_search%len_req
!
#ifdef PRISM_ASSERTION
      do i = 1, n_found
         if (index_found (i) <= 0 .and. &
             index_found (i) /= masked_out) exit
      end do
!
      if (i <= n_found) then
         print *, 'i, n_found, index_found', i, n_found, index_found (i)
         call psmile_assert ( __FILE__, __LINE__, &
                 'invalid index in found list')
      endif

      ip = 0
      do ipart = 1, search%npart
         do i = 1, len_nsend (ipart)
            if (indices_returned (ip+i) < 1 .or. &
                indices_returned (ip+i) > len_req(ipart)) exit
         end do
!
         if (i <= len_nsend (ipart)) then
            print *, 'ipart, ip, i, nreq, indices_returned(ip+i)', &
                     ipart, ip, i, len_req(ipart),                 &
                     indices_returned (ip+i)
            call psmile_assert ( __FILE__, __LINE__, &
                     'invalid index in indices_returned')
         endif
!
         ip = ip + len_nsend (ipart)
      end do
#endif
!
#ifdef DEBUGX
      do i = 1, n_found
         print *, 'i index_found', i, index_found(i)
      end do

      code = 1
      k = 0
      do i = 1, n_corners

         do n = 1, n_send
            if (IAND (found(n), code) /= 0) then
               k = k + 1
               print *, 'found: n', n, code, index_found (k)
            endif
         end do
         code = code * 2
      end do ! i

      ip = 0
      do ipart = 1, search%npart
         do i = 1, len_nsend (ipart)
               print *, 'ip i indices_returned', ip, i, indices_returned (ip+i), found (ip+i)
! indices_req(indices_returned (ip+i))
         end do
            ip = ip + len_nsend (ipart)
      end do
#endif
!
!===> Allocate array in order to store the indices of
!     of the compact list of "interpolation bases found"
!     (cf. PSMILe_Trili_3d_extra_off)
!
      Allocate (corner2ind (n_send, n_corners), stat = ierror)
      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = n_send * n_corners

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'corner2ind', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      corner2ind (:, :) = undef
!
!===> Decode mask "found" of points found
!     (cf. PSMILe_Trili_3d_extra_off)
!
      code = 1
      k = 0
      do i = 1, n_corners

         do n = 1, n_send
            if (IAND (found(n), code) /= 0) then
!
!===> ...... Required "i-th" interpolation base was found.
!            Corresponding value is stored in position "index_found(k)".
!            index_found (k) == masked_out : Point is masked out
!
               k = k + 1
               corner2ind (n, i) = index_found (k)
!
            endif
         end do
!
!===> ... next interpolation base
!
         code = code * 2
      end do ! i
!
#ifdef PRISM_ASSERTION
      if (k /= n_found) then
         print *, 'k, n_found', k, n_found
         call psmile_assert ( __FILE__, __LINE__, &
                 'k /= n_found')
      endif
#endif
!
!===> Remove points from code "found" which are masked out
!
      if (use_how /= PSMILE_novalue) then
!
         code = 1
         do i = 1, n_corners
            do n = 1, n_send
               if (corner2ind (n,i) == masked_out) then
                  found (n) = found (n) - code
               endif
            end do
            code = code * 2
         end do ! i
      endif
!
!===> For each part
!
!     indices_req = Indices of the extra points in entire (!) list of
!                   all points to be searched.
!     required    = Code for the extra points
!                   which are additionally required.
!
      ip = 0
      do ipart = 1, search%npart
         if (len_req (ipart) > 0 .and. &
             len_nsend (ipart) > 0) then
!
            indices_req => extra_search%indices_req(ipart)%vector
            required    => extra_search%required   (ipart)%vector
!
!===> ... Apply rules for points masked out
!
            if (use_how == PSMILE_novalue) then
!
!                 If one neighbor is invalid reject all 
!                 by moving all neighbors out of the valid range.

               do i = 1, n_corners
!cdir vector
                  do n = 1, len_nsend (ipart)
                     if (corner2ind (ip+n,i) == masked_out) then
!                       This point is not required any more
                        required (indices_returned (ip+n)) = 0
                        neighbors_3d (1, indices_req(indices_returned (ip+n)), &
                                      1:n_corners) = outside
                     endif
                  end do
               end do
            else 
!
! HR: This should be work also for "use_how == PSMILE_novalue"
!     "psmile_neigh_extra_points" should do the correct things in this case.
!     TODO: TEST whether this statement is really correct.
!
               do i = 1, n_corners
!cdir vector
                  do n = 1, len_nsend (ipart)
#ifdef DEBUG
                        print *, 'indices_returned', ip+n, indices_returned (ip+n), indices_req(indices_returned (ip+n))
                        print *, 'corner info', i, corner2ind (ip+n,i)
#endif
                     if (corner2ind (ip+n,i) == masked_out) then
                        neighbors_3d (1, indices_req(indices_returned (ip+n)), &
                                      i) = outside
#ifdef DEBUG_TRACE
                        if (indices_req(indices_returned (ip+n)) == ictl) then
                           print *, '#### psmile_add_points_found outside: ', outside
                           print *, '#### psmile_add_points_found ictl: set outside for ictl', ictl
                           print *, '#### psmile_add_points_found corner ', i
                           print *, ' -> neigh 1 ', neighbors_3d (1,indices_req(indices_returned (ip+n)),:)
                           print *, ' -> neigh 2 ', neighbors_3d (2,indices_req(indices_returned (ip+n)),:)
                           print *, ' -> neigh 3 ', neighbors_3d (3,indices_req(indices_returned (ip+n)),:)
                        endif
#endif
                     endif
                  end do
               end do
            endif
!
!===> ... Get common points required by "found" and "required" and
!         remove them from the coded list.
!
!cdir vector
            do n = 1, len_nsend (ipart)
               iband (n) = IAND (found (ip+n), required(indices_returned (ip+n)))
            end do
!
!cdir vector
            do n = 1, len_nsend (ipart)
               required(indices_returned (ip+n)) = &
               required(indices_returned (ip+n)) - iband (n)
            end do
!
!===> ... Set
!         neighbors_3d (1, i) = "extra_search%global_marker" in order to mark
!         that this is a point out of list of globally searched bases.
!         neighbors_3d (2, i) = Index in list of globally searched bases.
!
            code = 1
            do i = 1, n_corners
!
!cdir vector
               do n = 1, len_nsend (ipart)
                  if (IAND (iband(n), code) /= 0) then
#ifdef DEBUG_TRACE
                     if (indices_req(indices_returned (ip+n)) == ictl) then
                        print *, '#### ictl:', ictl, ' set global marker and', corner2ind(ip+n, i)
                        print *, 'i, ip, indices_returned', &
                                  i, ip, indices_returned (ip+n), iband (n)
                     endif
#endif
                     neighbors_3d (1,                              & 
                        indices_req(indices_returned (ip+n)), i) = &
                        extra_search%global_marker
!
                     neighbors_3d (2,                              & 
                        indices_req(indices_returned (ip+n)), i) = &
                        corner2ind(ip+n, i)
                  else
#ifdef DEBUG_TRACE
                     if (indices_req(indices_returned (ip+n)) == ictl) then
                        print *, '#### ictl:', ictl, ' global marker not set'
                        print *, 'i, ip, indices_returned', &
                                  i, ip, indices_returned (ip+n), iband (n), code
                     endif
#endif
                  endif
               end do
!
!===> ... next interpolation base
!
               code = code * 2
            end do
!
#ifdef PRISM_ASSERTION
assert1:    do i = 1, n_corners
!cdir vector
               do n = 1, len_nsend (ipart)
                  if (neighbors_3d (1, indices_req(indices_returned (ip+n)), i) &
                      == extra_search%global_marker .and. &
                      neighbors_3d (2, indices_req(indices_returned (ip+n)), i) &
                      <= 0)  exit assert1
               end do
            end do assert1
!
            if (i <= n_corners) then
               print *, 'ipart, i, n', ipart, i, n, corner2ind(ip+n, i), &
                                       found(ip+n), iband(n)
               call psmile_assert ( __FILE__, __LINE__, &
                        'invalid index for interpolation base generated')
            endif
#endif
         endif
!
         ip = ip + len_nsend (ipart)
      end do ! ipart 
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

9990 format (1x, a, ': psmile_add_points_found: n_send ', i6, &
             '; n_found ', i6)
9980 format (1x, a, ': psmile_add_points_found: eof ierror =', i3)

#endif /* VERBOSE */

#ifdef DEBUG
#endif

      end subroutine PSMILe_Add_points_found