psmile_neigh_global_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_global_points
!
! !INTERFACE:

      subroutine psmile_neigh_global_points (search, extra_search,     &
                      mask_array, mask_shape, mask_available, use_how, &
                      grid_id, grid_valid_shape,                       &
                      neighbors_3d, nloc, n_corners, len_cpl, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_global_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

!     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_id
!
!     Grid identifier of the source grid
!     (on which the target points are searched)
!
      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%search_data%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_global_points;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!     ... for grids
!
      Integer                         :: grid_type
!
!     ... for locations searched
!
      Integer                         :: i, j, n, nbeg, nend, nprev
!
!     ... for locations controlled
!
      Integer                         :: code, nreq
      Integer                         :: outside

!
      Integer, Pointer                :: required (:)
      Integer, Pointer                :: indices_req (:)
      Integer, Allocatable            :: mask_cell (:)
      Integer, Allocatable            :: search_required (:)
!
!     ... for partitions
!
      Integer                         :: ipart

!     ... for error handling
!
      Integer, parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_global_points" determines the points (indices)
! which require global 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_global_points.F90 3174 2011-05-04 09:41:43Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_global_points.F90 3174 2011-05-04 09:41:43Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
!===> Initialization
!
      ierror  = 0
!
#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
#endif
!
      if (n_corners > 31) then
         ierror = PRISM_Error_internal
         ierrp (1) = n_corners

         call psmile_error ( ierror, &
                            'Number of corrners too large', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      endif
!
!------------------------------------------------------------------------
!     Look for points required in global search
!     Value of "search_required" contains a code (2*(n-1)) for the
!     interpolation points to be searched in the global search.
!
!     search_required (i) = 0: No search required for i-th point
!------------------------------------------------------------------------
!
! NOTE: In the case of hyperplanes (length(j) == 1), 
!       interpolation bases may be multiple requested since
!       neighbors (:, i, n) may be the identical for 
!       different corners.
!       TODO: Remove multiple entries. Should be possibly done
!             before sending the data to adjoining block.
!
      Allocate (search_required (nloc), STAT = ierror)

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

         ierror = PRISM_Error_Alloc

         call psmile_error ( ierror, 'search_required', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      search_required (:) = 0
!
!===> Code points for which global searched is required
!
      grid_type = Grids(grid_id)%grid_type
!
      if (grid_type == PRISM_Reglonlatvrt     .or. &
          grid_type == PRISM_Irrlonlat_Regvrt .or. &
          grid_type == PRISM_Irrlonlatvrt     .or. &
          grid_type == PRISM_Gaussreduced_Regvrt) then
!
         code = 1
            do n = 1, n_corners
!cdir vector
               do i = 1, nloc
               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
                  search_required (i) = search_required (i) + code
               end if
               end do
!
!===> ... next interpolation base (corner)
!
            code = code * 2
            end do
!
!-----------------------------------------------------------------------
!     Control Mask values
!
! 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.
!-----------------------------------------------------------------------
!
         if (mask_available .and. use_how == PSMILE_novalue) then
!
!===> ... Allocate and initialize work vector "mask_cell" with
!         mask_cell (i) = Number of neighbour points of cell "i"
!                         which are masked out.
!
            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
!
            mask_cell (:) = 0

!           Get corners which are flagged as false

            code = 1

            do n = 1, n_corners
#ifdef DEBUG_TRACE
               if (ictl > 0 .and. ictl <= nloc) then
                  if (neighbors_3d(1, ictl, n) >= grid_valid_shape(1,1) .and. &
                      neighbors_3d(1, ictl, n) <= grid_valid_shape(2,1) .and. &
                      neighbors_3d(2, ictl, n) >= grid_valid_shape(1,2) .and. &
                      neighbors_3d(2, ictl, n) <= grid_valid_shape(2,2) .and. &
                      neighbors_3d(3, ictl, n) >= grid_valid_shape(1,3) .and. &
                      neighbors_3d(3, ictl, n) <= grid_valid_shape(2,3) ) then
                     print *, 'DEBUG_TRACE: psmile_neigh_global_points: ictl',   &
                              ictl, n, neighbors_3d(:, ictl, n),    &
                              mask_array (neighbors_3d(1, ictl, n), &
                                          neighbors_3d(2, ictl, n), &
                                          neighbors_3d(3, ictl, n))
                 endif
              endif
#endif

!cdir vector
               do i = 1, nloc
               if (IAND (search_required (i), code) == 0) then
!
!                 neighbors_3d (:, i, n) is within valid grid shape
!                 See above
!
                  if ( .not. mask_array (neighbors_3d(1, i, n), &
                                         neighbors_3d(2, i, n), &
                                         neighbors_3d(3, i, n)) ) &
                     mask_cell(i) = mask_cell (i) + 1
               endif
               end do
!
!===> ... next interpolation base (corner)
!
            code = code * 2
            end do ! n
!
!     Note: "outside" must be greater than
!           "grid_valid_shape (2,1) + maximal with of interpolation stencil"
!           in order to enable global search.
!               
          if ((grid_type == PRISM_Gaussreduced_Regvrt) .and. &
              associated(grids(grid_id)%reduced_gauss_data%remote_index)) then
            outside = grid_valid_shape(2,1) + &
                      size(Grids(grid_id)%reduced_gauss_data%remote_index) + 11
          else
            outside = grid_valid_shape(2,1) + 11
          endif
!
!===> How should masked values be handled ?
!
!       if ( use_how == 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)
!
               do n = 1, n_corners
!cdir vector
                  do i = 1, nloc
                  if ( mask_cell(i) > 0 ) then
                     neighbors_3d(1, i, n) = outside
                  endif
                  enddo
               enddo
!
!===>          Remove all points from global search for which
!              already a masked out value was found.
!
               do i = 1, nloc
                  if ( mask_cell(i) > 0) then
                     search_required (i) = 0
                  endif
               end do
!
!           endif ! if (use_how == PSMILE_novalue)
!
            Deallocate (mask_cell)

         endif

      else 

         ierrp (1) = grid_type
         ierror = PRISM_Error_Internal

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

      endif
!
!===> Store cells/points which require global search
!     Note: The index in the entire (!) vector "neighbors_3d" is stored
!           and not the relative index in the partition "ipart"
!     ??? Lohnt es sich das fuer jede Partition abzuspeichern ?
!     ??? Geht das nicht fuer alle Punkte ?
!
      nprev = 0
!
      do ipart = 1, search%search_data%npart
         if (len_cpl (ipart) > 0) then
            nbeg = nprev + 1
            nend = nprev + len_cpl (ipart)
!
            nreq = 0
!cdir vector
            do i = nbeg, nend
               if (search_required(i) /= 0) nreq = nreq + 1
            end do

            if (nreq > 0) then
               Allocate (indices_req (nreq), required (nreq), stat = ierror)
               if (ierror > 0) then
                  ierrp (1) = ierror
                  ierrp (2) = nreq * 2

                  ierror = PRISM_Error_Alloc

                  call psmile_error ( ierror, 'indices_req, required', &
                                      ierrp, 2, __FILE__, __LINE__ )
                  return
               endif
!
               j = 0
!cdir vector
               do i = nbeg, nend
                  if (search_required(i) /= 0) then
                     j = j + 1
                     indices_req (j) = i
                     required (j) = search_required (i)
                  endif
               end do
!
!===> ...... Save global search data in structure "extra_search"
!
               extra_search%indices_req(ipart)%vector => indices_req
               extra_search%required   (ipart)%vector => required
               extra_search%len_req    (ipart)        = nreq
!
               extra_search%n_req = extra_search%n_req + nreq

#ifdef DEBUG_TRACE
               do j = 1, nreq
                  if (ictl == indices_req (j)) exit
               end do
!
               if (j <= nreq) then
                  print *, 'DEBUG_TRACE: psmile_neigh_global_points: ipart, ictl', &
                             ipart, ictl
                  print *, 'DEBUG_TRACE: index in indices_req :', j, &
                              ', required:', required(j)
!
                  do n = 1, n_corners
                     print *, 'DEBUG_TRACE: ng', neighbors_3d(:, ictl, n)
                  end do
               endif
#endif 
!
#ifdef DEBUGX
               print *, 'ipart, nreq ', ipart, nreq
               do j = 1, nreq
                  i = indices_req (j)
                  print *, 'j, i, required', j, i, required(j)
!
                  do n = 1, n_corners
                     print *, 'ng', neighbors_3d(:, i, n)
                  end do
               end do
#endif
            endif ! nreq > 0
!
!===> ... next partition
!
         nprev = nend
         endif ! len_cpl  > 0
      end do ! ipart
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), extra_search%n_req

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990 format (1x, a, ': psmile_neigh_global_points:')
9980 format (1x, a, ': psmile_neigh_global_points: eof, nreq', i6)

#ifdef PRISM_ASSERTION
#endif

      end subroutine psmile_neigh_global_points