psmile_trili_gauss2_extra.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! Copyright 2011, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: psmile_trili_gauss2_extra
!
! !INTERFACE:

      subroutine psmile_trili_gauss2_extra (search, grid_id, mask_array, &
                        mask_shape, mask_available, ibuf, len_item,      &
                        n_send, num_neigh, ierror)
!
! !USES:
!
      use PRISM
      use PSMILe, dummy_interface => psmile_trili_gauss2_extra
      use psmile_grid_reduced_gauss

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer,                      Intent (In)    :: grid_id

      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)     :: len_item
!
!     Number of pieces of information in integer buffer "ibuf".
!        1:3: (i,j,k) index of lower left corner
!        4  : Index in
!        5  : Code for points required
!
      Integer,                     Intent (In)    :: n_send
!
!     Number of extra points sent

      Integer,                     Intent (In)    :: num_neigh
!
!     (Maximal) Number of interpolations bases required..
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Enddef_global_search), Intent (InOut) :: search

!     Data on the points to be searched

      Integer,                     Intent (InOut) :: ibuf (len_item, n_send)
!
!     Integer buffer containing the packed integer data.
!     On return, ibuf (5, *) contains the code of the points found.
!
! !OUTPUT PARAMETERS:
!
      Integer,                     Intent (Out)   :: ierror

!     Returns the error code of PSMILe_Search_donor_extra;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS:
!
! n_corners_3d = Number of interpolation bases in 3d case
!
! masked_out   = Hash value for a point which is masked out
!
      Integer, Parameter           :: n_corners_3d = 2 * 2 * 2
!
      Integer, Parameter           :: masked_out = 0
!
! !LOCAL VARIABLES
!
!     ... loop variables
!
      Integer                      :: j, n, i, k
!
!     ... For Grids
!
      Type (Grid), Pointer         :: grid_info
      Integer,     Pointer         :: grid_valid_shape (:, :)
!
!     ... For Gaussreduced Grids
!
      Logical                      :: virtual_cell_available
!
!     ... For searched points of control volume
!
      Integer                      :: code, n_corners, outside
!
      Integer, Allocatable         :: locs (:, :)
      Integer, Allocatable         :: neighbours_3d_extra (:, :, :)
      Integer, Allocatable         :: hash (:)
!
!
!
      Logical                      :: missing_points(num_neigh, n_send)
      Integer                      :: decode_mask(num_neigh)
      Integer                      :: num_missing_points

      integer, parameter           :: shift_3d(ndim_3d, 9) = 
                                          reshape ((/-1,-1,0, 0,-1,0, +1,-1,0,
                                                     -1, 0,0, 0, 0,0, +1, 0,0,
                                                     -1,+1,0, 0,+1,0, +1,+1,0/), 
                                                   (/ndim_3d,9/))
!
!     ... for error parameters
!
      Integer, parameter           :: nerrp = 1
      Integer                      :: ierrp (nerrp)
!
#ifdef DEBUG_TRACE
      Integer                      :: ictl_req, ind_loc (n_corners_3d)
#endif
!
! !DESCRIPTION:
!
! Subroutine "psmile_trili_gauss2_extra" performs the
! additional (global) search for grids of type "PRISM_Gaussreduced_regvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.05    H. Ritzdorf  created as psmile_trili_gauss2_extra_dble
! 17.02.11    M. Hanke     simplified code and corrected shifting, removed difference
!                          between dble and real implementation
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_trili_gauss2_extra.F90 3023 2011-03-17 10:21:20Z hanke $
! $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_trili_gauss2_extra.F90 3023 2011-03-17 10:21:20Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), grids(grid_id)%comp_id

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
      n_corners = num_neigh
!
!===> Get grid info
!
!  grid_valid_shape (2, ndim_3d) : Specifies the valid block shape
!                                  for the "ndim_3d"-dimensional block
!                                  (without halo/overlap regions)
!
      grid_info => Grids(grid_id)
!
      grid_valid_shape => grid_info%grid_shape
!
      virtual_cell_available = len_item > ndim_3d + 2 ! unschoen
!
      if (grid_info%grid_type /= PRISM_Gaussreduced_Regvrt) then
         print *, "grid_type", grid_info%grid_type
         call psmile_assert (__FILE__, __LINE__, &
                   "This routine is not designed for this grid type")
      endif
!
      outside = grid_valid_shape (1, 1) - 11
!
#ifdef PRISM_ASSERTION

      if (num_neigh < 1 .or. num_neigh > n_corners_3d) then
         print *, 'num_neigh', num_neigh, n_corners_3d
         call psmile_assert ( __FILE__, __LINE__, &
                'Invalid value for num_neigh (number of interpolation bases')
      endif
!
!===> Control input values contained in "ibuf"
!
!cdir vector
      do n = 1, n_send
         if (ibuf (2, n) < 1 .or. ibuf (2, n) > grid_info%nbr_latitudes) exit
      end do
!
      if (n <= n_send) then
         print *, "i,j,k, gnbr_lats", ibuf(1:3,n), grid_info%nbr_latitudes
         call psmile_assert (__FILE__, __LINE__, &
                   "second (virtual) index out of range")
      endif
!
!cdir vector
      do n = 1, n_send
         if (ibuf (1, n) < 1 .or. &
             ibuf (1, n) > grid_info%nbr_points_per_lat(ibuf(2,n))) exit
      end do
!
      if (n <= n_send) then
         print *, "i,j,k, points_per_lat", &
                  ibuf(1:3,n), grid_info%nbr_points_per_lat(ibuf(2,n))
         call psmile_assert (__FILE__, __LINE__, &
                   "first (virtual) index out of range")
      endif
#endif

#ifdef DEBUG_TRACE
!cdir vector
      do ictl_req = 1, n_send
         if (ibuf (4, ictl_req) == 1) exit
      end do
#endif

      ! 1st we rebuilt the neighbours_3d (global 3d representation) array
      ! for the source locations received

      allocate (neighbours_3d_extra (ndim_3d, n_send, 4), stat=ierror)
      if (ierror /= 0) then
         call psmile_error ( PRISM_Error_Alloc, 'locs', (/ndim_3d * n_send * 4/), 1, &
                 __FILE__, __LINE__ )
         return
      endif

      if (virtual_cell_available) then
         ! for all locations
         do n = 1, n_send

            if (ibuf(6, n) == 0) then

               neighbours_3d_extra(:,n,1:4) = &
                  psmile_gauss_get_bili_stencil (grid_id, ibuf(1:3,n))

            else

               neighbours_3d_extra(:,n,1:4) = &
                  psmile_gauss_shift_bili_stencil (grid_id, ibuf(1:3,n), shift_3d(:,ibuf(6, n)))

            endif
         enddo ! n
      else
         neighbours_3d_extra(:,n,1:4) = psmile_gauss_get_bili_stencil (grid_id, ibuf(1:3,n))
      endif
!
!  decode which points are missing
!
      do i = 1, num_neigh
         decode_mask(i) = ishft (1, i - 1)
      enddo

      do n = 1, n_send
         missing_points(:, n) = iand (decode_mask(:), ibuf(5, n)) /= 0
      enddo
!
!  how many points are missing?
!
      num_missing_points = count (missing_points)
#ifdef PRISM_ASSERTION
      if (num_missing_points <= 0) then
         call psmile_assert ( __FILE__, __LINE__, &
                             'Number of points to be searched == 0 ?!?')
      endif
#endif

#ifdef DEBUGX
      print *, 'code for locations to be searched (n_send =', num_missing_points, ')'
      do n = 1, n_send
         print *, n, 'srcloc =', ibuf(1:3,n), 'code =', ibuf (5, n)
         print *, '   decoded     =', missing_points(:,n)
      end do
#endif
!
#ifdef DEBUG_TRACE
      if (ictl_req <= n_send) then
         print *, 'code for indices_req', ibuf (4, ictl_req)
         print *, 'ictl_req', ictl_req, 'srcloc', ibuf(1:3,ictl_req), &
                  'code', ibuf (5, ictl_req), missing_points(:, ictl_req)
      endif
#endif
!
!  Allocate memory for the indices of the missing points
!
      Allocate (locs (ndim_3d, num_missing_points), stat = ierror)
      if (ierror /= 0) then
         call psmile_error ( PRISM_Error_Alloc, 'locs', (/ndim_3d * num_missing_points/), 1, &
                 __FILE__, __LINE__ )
         return
      endif

      k = 0

      do i = 1, n_corners
!cdir vector
         do n = 1, n_send
            if (missing_points (i,n)) then

               k = k + 1

               locs (:, k) = neighbours_3d_extra(:, n, iand (i-1, 3)+1)

               if (grid_valid_shape(2,3) - grid_valid_shape(1,3) >= ishft(n_corners, -2) - 1 .and. &
                   n_corners > 4) then
                  locs (3, k) = locs (3, k) + ishft (i-1, -2)
               endif
            end if
         end do ! n
      end do ! i
!
#ifdef PRISM_ASSERTION
      if (k /= num_missing_points) then
         print *, 'k, num_missing_points', k, num_missing_points
         call psmile_assert ( __FILE__, __LINE__, &
                'Inconsistent values of k and num_missing_points; must be same value')
      endif
#endif

#ifdef DEBUG_TRACE
      if (ictl_req <= n_send) then
         ind_loc (:) = 0
         n = 0
         do j = 1, n_corners
!cdir vector
            do i = 1, n_send
               if (missing_points (j,i)) then
                  n = n + 1
                  if (i == ictl_req) ind_loc (j) = n
               end if
            end do ! i
         end do ! j
!
         print *, 'locs searched for indices_req', ibuf (4, ictl_req)
         do j = 1, n_corners
            n = ind_loc (j)
            if (n > 0) &
            print *, 'ictl_req', ictl_req, ', n ', n, &
                           ': locs', locs (1:ndim_3d, n)
         end do
      endif
#endif

!
!===> Transform indices from global (virtual) 3d indices
!                       into local  (virtual) 3d indices
!
      locs(1,:) = psmile_gauss_3d_to_local_1d (grid_id, locs(:,:), &
                                               num_missing_points, outside)
      locs(2,:) = 1
      !locs(3,:) = locs(3,:)
!
!===> Create compact list of points to be sent and
!     return hash value "hash" (index in fortran order) for each location
!
      Allocate (hash (num_missing_points), stat = ierror)
      if (ierror /= 0) then
         ierrp (1) = num_missing_points

         ierror = PRISM_Error_Alloc

         call psmile_error ( ierror, 'hash', ierrp, 1, &
                 __FILE__, __LINE__ )
         return
      endif
!
      call psmile_hash_extra (search, locs, hash, num_missing_points,                &
              mask_array, mask_shape, mask_available, grid_valid_shape, &
              ierror)
      if (ierror > 0) return
!
#ifdef DEBUGX
      print *, 'hash values and locs ', grid_valid_shape
      n = 0
      do j = 1, n_corners
!cdir vector
         do i = 1, n_send
            if (missing_points (j,i)) then
               n = n + 1
               print *, 'i, j, n, hash(n), locs(:,n) ', &
                         i, j, n, hash(n), locs(:,n)
            end if
         end do ! i
      end do ! j
#endif
!
!===> Save coded mask of points found in "ibuf"
!     ??? Ob das so clever ist ?
!
      if (search%n_liste > 0) then
         ibuf (5, 1:n_send) = 0
!
         n = 0
         code = 1
         do j = 1, n_corners
!cdir vector
            do i = 1, n_send
               if (missing_points (j,i)) then

                  n = n + 1
!
                  if (hash (n) >= 0) then
!
!===> ...... Required "j-th" interpolation base was found
!            (may be masked out (hash(n) == masked_out).
!            Corresponding value is stored in position "index_found(n)".
!
                     ibuf (5, i) = ibuf (5, i) + code
                  endif
!
               endif
            end do ! i

            code = code * 2
         end do ! j

#ifdef PRISM_ASSERTION
         if (n /= num_missing_points) then
            print *, 'num_missing_points, n', num_missing_points, n
            call psmile_assert ( __FILE__, __LINE__, &
                                'n /= num_missing_points')
         endif
#endif
!
#ifdef DEBUGX
         print *, 'code for locations found', n_send
         do i = 1, n_send
            print *, 'code', i, ibuf (5, i), missing_points(:, i)
         end do
#endif
!
#ifdef DEBUG_TRACE
         if (ictl_req <= n_send) then
            print *, 'code for locations found for indices_req', ibuf (4, ictl_req)
            print *, 'ictl_req', ictl_req, ': code', ibuf (5, ictl_req), missing_points(:, ictl_req)
         endif
#endif

      endif
!
!===> Deallocate memory allocated
!
      Deallocate (locs, hash, neighbours_3d_extra)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror, search%n_found, n_send

      call psmile_flushstd
#endif /* VERBOSE */
!
      return
!
!  Formats:
!
#ifdef VERBOSE
9990 format (1x, a, ': psmile_trili_gauss2_extra: comp_id =', i3)
9980 format (1x, a, ': psmile_trili_gauss2_extra: eof', &
                    ', ierror =', i4, ', n_found =', i8, i8)
#endif /* VERBOSE */

      end subroutine psmile_trili_gauss2_extra