psmile_neigh_extra_points.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_Neigh_extra_points
!
! !INTERFACE:

      subroutine psmile_neigh_extra_points (search, extra_search, &
               mask_array, mask_shape, mask_available, use_how,   &
               grid_valid_shape,                                  &
               neighbors_3d, nloc, n_corners, len_cpl, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_extra_points
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Type (Enddef_search), Intent (In) :: search

!     Info's on coordinates to be searched
!
      Integer,         Intent (In)    :: mask_shape (2, ndim_3d)

!     Dimension of (method) mask array

      Logical,         Intent (In)    :: mask_array (mask_shape (1,1): 
                                                     mask_shape (2,1), 
                                                     mask_shape (1,2): 
                                                     mask_shape (2,2), 
                                                     mask_shape (1,3): 
                                                     mask_shape (2,3))
!     Mask of the method

      Logical,         Intent (In)    :: mask_available

!     Is mask specified in array "mask_array" ?
!     mask_available = .false. : Mask is not available
!                      .true.  : Mask is     available

      Integer,         Intent (In)    :: use_how(3)

!     Information about how to apply the mask
!        PSMILE_novalue   : If one neighbor is unvalid reject all.
!        PSMILE_tneighbour: Reject all unvalid neighbors but keep the valid.
!        PSMILE_nneighbour: == PSMILE_undef

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

!     Specifies the valid block shape for the "ndim_3d"-dimensional block
!
      Integer,         Intent (In)    :: nloc
!
!     Total number of locations to be transferred.
!     Number of locations for which extra points should be searched.
!     Second dimension of neighbors array "neighbors_3d"
!
      Integer,         Intent (In)    :: n_corners
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of interpolation bases.
!
      Integer, Intent (In)            :: len_cpl (search%npart)

!     Number of points to be sent to the coupler for each partition.
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Extra_search_info), Intent (InOut) :: extra_search
!
!     Number of locations where
!       (*) required mask values were not true
!
      Integer,         Intent (InOut) :: neighbors_3d (ndim_3d, nloc, n_corners)

!     Indices of neighbor locations (interpolation bases)
!
! !OUTPUT PARAMETERS:

      Integer,         Intent (Out)   :: ierror

!     Returns the error code of PSMILE_Neigh_extra_points;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: grid_id
      Integer, Pointer                :: recv_list(:)
      Integer                         :: search_mode
      Integer                         :: i, j, n, nbeg, nend
!
!     ... for locations controlled
!
      Integer                         :: n_extra, outside
!
      Integer, Pointer                :: indices (:)
      Logical, Allocatable            :: mask_cell (:)
!
!     ... for partitions
!
      Integer                         :: ipart
!
!     ... for error handling
!
      Integer, parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_extra_points" determines the points (indices)
! which require extra search (nearest neighbour search).
! The extra points are searched the neighbour indices (interpolation
! bases) specified in
!
!      neighbour (1:ndim_3d, 1:nloc, n_corners).
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 15.03.06    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_extra_points.F90 2991 2011-02-25 08:34:51Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_extra_points.F90 2991 2011-02-25 08:34:51Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
!===> Initialization
!
      ierror  = 0
      search_mode = use_how(1)
      grid_id = use_how(3)
!
!     Note: "outside" must be greater than
!           "grid_valid_shape (2,1) + maximal with of interpolation stencil"
!            in order to enable global search.
!
      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
!
#ifdef PRISM_ASSERTION
!
      if (nloc /= Sum(len_cpl)) then
         print *, "nloc, Sum(len_cpl)", nloc, Sum (len_cpl)
         call psmile_assert (__FILE__, &
              __LINE__, "ncpl == SUM(len_cpl) expected")
      endif

      if (extra_search%n_extra /= 0) then
         call psmile_assert (__FILE__, &
              __LINE__, "extra_search%n_extra /= 0")
      endif
#endif
!
!-----------------------------------------------------------------------
!     Control Mask values
!
!     The searched point is located in donor cell
!     (srcloc(:,i), srcloc(:,i)+1).
!-----------------------------------------------------------------------
!
! PSMILE_novalue  : If one neighbor is unvalid reject all.
!                   This will cause the transformer to perform
!                   bi/trilinear interpolation for those points
!                   only for which all neighbors are valid.
!
! PSMILE_tneighbour: Reject all unvalid neighbors but keep the valid.
!                    This will cause the transformer to perform
!                    bi/trilinear interpolation for those points
!                    only for which all neighbors are valid.
!                    Those for which less than "n_corners" valid neighbors
!                    have been found will be calculated using weighted
!                    distances.
!
! PSMILE_nneighbour: == PSMILE_undef [default]
!                    Bi/trilinear interpolation for those points for which
!                    all "n_corners" points are valid.
!                    For those for which less than "n_corners" valid
!                    neighbors were found those neighbors are kept.
!                    For the "n_extra" target points for which 0 valid
!                    neighbors were found we invoke the search for
!                    one single nearest neighbor.
!
! NOTE: If a interpolation base is outside of the physical boundary,
!       the base points is handled such as the mask would be false.
!
!-----------------------------------------------------------------------
!
!     set search mode PSMILE_nneighbour to default
!
#ifdef DEBUG
8990  format (1x, a, ': psmile_neigh_extra_points: search_mode', i8, &
              '; mask_available ', l1)
      print 8990, trim(ch_id), search_mode, mask_available
#endif
!
      if ( search_mode == PSMILE_nneighbour ) search_mode = PSMILe_Undef

#ifdef DEBUG_TRACE
         if ( ictl > 0 .and. ictl <= nloc ) then
            print *, ' ictl neigh 1', neighbors_3d(1,ictl,:)
            print *, ' ictl neigh 2', neighbors_3d(2,ictl,:)
            print *, ' ictl neigh 3', neighbors_3d(3,ictl,:)
         endif
#endif
 if (mask_available) then

!        Ignore corners of (method) cells flagged as false by moving them
!        out of the valid range (see psmile_compact_neighbors_3d)

         do n = 1, n_corners
!cdir vector
            do i = 1, nloc

              if (neighbors_3d(1, i, n) /= extra_search%global_marker) then
!
!             search_mode == vneighbour, not implemented !
!
!             If the point was not found in regular domain,
!             mark it by "outside".
!
                  if ( neighbors_3d(1, i, n) < grid_valid_shape(1,1) .or. &
                       neighbors_3d(1, i, n) > grid_valid_shape(2,1) .or. &
                       neighbors_3d(2, i, n) < grid_valid_shape(1,2) .or. &
                       neighbors_3d(2, i, n) > grid_valid_shape(2,2) .or. &
                       neighbors_3d(3, i, n) < grid_valid_shape(1,3) .or. &
                       neighbors_3d(3, i, n) > grid_valid_shape(2,3) ) then
!
                       neighbors_3d(1, i, n) = outside
!
!             This corresponds to "search_mode == tneighbour".
!             Only the actual neighour point is moved outside
!             since the mask is .false..
!
                  else if ( neighbors_3d(1, i, n) >= mask_shape(1,1) .and. &
                            neighbors_3d(1, i, n) <= mask_shape(2,1) .and. &
                            neighbors_3d(2, i, n) >= mask_shape(1,2) .and. &
                            neighbors_3d(2, i, n) <= mask_shape(2,2) .and. &
                            neighbors_3d(3, i, n) >= mask_shape(1,3) .and. &
                            neighbors_3d(3, i, n) <= mask_shape(2,3) ) then
#ifdef DEBUG_TRACE
                     if ( i == ictl ) then
                          print *, ' ictl neigh & mask ', i, n, neighbors_3d(:, i, n), &
                                                    mask_array (neighbors_3d(1, i, n), &
                                                                neighbors_3d(2, i, n), &
                                                                neighbors_3d(3, i, n))
                     endif
#endif
                     if ( .not. mask_array (neighbors_3d(1, i, n),   &
                                            neighbors_3d(2, i, n),   &
                                            neighbors_3d(3, i, n)) ) &
                       neighbors_3d(1, i, n) = outside
                  endif

               endif
            enddo ! i
         enddo ! n
!
!    search_mode == tneighbour is already treated above. Now we check whether
!    additional neighbours have to be marked as not valid (outside) for
!        search_mode == novalue
!        search_mode == nneighbour | undef [default]  
!
         if ( search_mode == PSMILE_novalue .or. &
              search_mode == PSMILe_undef ) then
!
!===> ... Allocate and initialize work vector "mask_cell" with
!         mask_cell (i) = .true. : Mask of first point is .false.
!                                  ==> First index is out of range
!
            Allocate (mask_cell (nloc), STAT = ierror)
            if (ierror > 0) then
               ierrp (1) = ierror
               ierrp (2) = nloc

               ierror = PRISM_Error_Alloc

               call psmile_error ( ierror, 'mask_cell', &
                    ierrp, 2, __FILE__, __LINE__ )
               return
            endif
!
!cdir vector
            do i = 1, nloc
               mask_cell (i) = neighbors_3d(1, i, 1) == outside
            end do
!
            if ( search_mode == PSMILE_novalue ) then
!
!           Ignore cells with one or more neighbors flagged as false 
!           by moving all neighbors out of the valid range
!           (see psmile_compact_neighbors_3d)
!
!           mask_cell (i) = .true. : A mask entry of (method) cell is .false.
!
               do n = 2, n_corners
!cdir vector
                  do i = 1, nloc
                     mask_cell (i) = mask_cell (i) .or. &
                          neighbors_3d(1, i, n) == outside
                  end do
               end do

               do n = 1, n_corners
                  !cdir vector
                  do i = 1, nloc
                     if ( mask_cell(i) ) &
                        neighbors_3d(1, i, n) = outside
                  enddo
               enddo
!
            else ! if ( search_mode == PSMILe_undef ) then
!
!               mask_cell (i) = .true. : All mask entries of (method) cell
!                                        are .false.
!
               do n = 2, n_corners
!cdir vector
                  do i = 1, nloc
                     mask_cell (i) = mask_cell (i) .and. &
                          neighbors_3d(1, i, n) == outside
                  end do
               end do
!
#ifdef TEST_EXTRA_SEARCH
               print *, '### psmile_neigh_extra_points: Extra-Search for testing enabled'
               mask_cell(:) = .true.
#endif
!
!===> Store cells/points which require extra search
!     Note: The index in the entire vector "neighbors_3d" is stored
!           and not the relative index in the partition "ipart"
!
!===> ... Get number of points to be searched extra
!
               n_extra = count (mask_cell (:))
!
               if (n_extra > 0) then
!
!===> ...... Allocate memory required to store the indices and
!            code for the points required.
!            Note: The indices are stored as index in entire (!) 
!                  array "neighbors"
!
                  Allocate (indices (n_extra), stat = ierror)
                  if (ierror > 0) then
                     ierrp (1) = ierror
                     ierrp (2) = n_extra

                     ierror = PRISM_Error_Alloc

                     call psmile_error ( ierror, 'extra_search%indices', &
                          ierrp, 2, __FILE__, __LINE__ )
                     return
endif
!
!===> ...... Generate indices
!
#ifdef USE_PACK
                  indices (1: n_extra) = &
                       PACK ((/ (i, i=1,nloc) /), mask_cell)
#else
                  j = 0
!cdir vector
                  do i = 1, nloc
                     if (mask_cell (i)) then
                        j = j + 1
                        indices (j) = i
                     endif
                  end do
#endif /* USE_PACK */
!
!===> ...... Generate length per partition
!
                  nend = 0
                  do ipart = 1, search%npart
                     if (len_cpl (ipart) > 0) then
                        nbeg = nend + 1
                        nend = nend + len_cpl (ipart)
!
                        extra_search%len_extra (ipart) = &
                             count (mask_cell(nbeg:nend))
!
                     endif ! len_cpl > 0
                  end do ! ipart
!
#ifdef DEBUGX
                  do i = 1, nloc
                     if (mask_cell (i)) then
                        print *, 'i', i, neighbors_3d(:, i, 1)
                        do n = 1, n_corners
                           print *, 'ng', neighbors_3d(:, i, n), &
                                mask_array (neighbors_3d(1, i, n), &
                                neighbors_3d(2, i, n), &
                                neighbors_3d(3, i, n))
                        end do
                     endif
                  end do
#endif
!
!===> ...... Save extra search data in structure "extra_search"
!
                  extra_search%indices => indices
!
                  extra_search%n_extra = n_extra

               endif ! (n_extra > 0)
            endif ! if (earch_mode == PSMILE_novalue)
!
            Deallocate (mask_cell)
            !
         endif ! if (search_mode == PSMILE_novalue .or. ...

      endif ! if (mask_available)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), extra_search%n_extra

!     print *, extra_search%len_extra (:)
      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990  format (1x, a, ': psmile_neigh_extra_points:')
9980  format (1x, a, ': psmile_neigh_extra_points: eof, n_extra', i6)

    end subroutine psmile_neigh_extra_points