psmile_neigh_tricu_irreg2.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_tricu_irreg2
!
! !INTERFACE:

      subroutine psmile_neigh_tricu_irreg2 (                 &
               grid_valid_shape, interpolation_mode, cyclic, &
               srcloc, srclocz, nlocs, nloc, npreva,         &
               neighbors_3d, num_neigh, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_tricu_irreg2

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: nloc
!
!     Total number of locations to be transferred
!
      Integer, Intent (In)            :: nlocs (ndim_2d)
!
!     Number of locations in each direction
!
      Integer, Intent (In)            :: srcloc  (ndim_2d, nlocs (1))
      Integer, Intent (In)            :: srclocz (nlocs (2))
!
!     Indices of the (method) grid cell in which the point was found.

      Integer, Intent (In)            :: npreva
!
!     Number of locations for which neighbours were already searched.
!
      Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)

!     Specifies the valid block shape for the "ndim_3d"-dimensional block
!
      Integer, Intent (In)            :: interpolation_mode
!
!     Specifies the interpolation mode with
!        interpolation_mode = 3 : Tri-linear or bilinear/linear interpolation
!                                 (determines 8 corners)
!        interpolation_mode = 2 : Bilinear interpolation
!                                 (determines 4 corners)
!        interpolation_mode = 1 : Linear   interpolation
!                                 (determines 2 corners)
!
      Logical, Intent (In)            :: cyclic (ndim_3d)
!    
!     Does the block have cyclic coordinates in the corresponding direction ?
!
      Integer, Intent (In)            :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d"
!
! !OUTPUT PARAMETERS:
!
!     Location of indices for neighbors_3d
!
!           13   14   15   16
!
!     ^      9   10   11   12
!     |
!     |      5    6    7    8
!
!     j      1    2    3    4
!
!       i ---->
!
      Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)

!     Indices of neighbor locations (interpolation bases)

      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_Neigh_tricu_irreg2;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  N_CORNERS_3D = Number of corners for a 3-d interpolation
!               = 2 ** ndim_3d
!  N_CORNERS_2D = Number of corners for a 2-d interpolation
!               = 2 ** ndim_2d
!
      Integer, parameter              :: n_corners_3d = 2 * 2 * 2 * 8
      Integer, parameter              :: n_corners_2d = 2 * 2 * 4
      Integer, parameter              :: n_corners_1d = 2 * 2
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: i, j, k, n
      Integer                         :: nbeg, nend, nprev
      Integer                         :: lenijk, n_corners
!
!     ... For locations controlled
!
      Integer                         :: length (ndim_3d)
      Integer                         :: nca (ndim_3d)
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_tricu_irreg2" searches the "n_corners"
! (8, 4 or 2) neighbours (3d-indices of the interpolation bases) for the
! selected interpolation type (tri-, bi-linear or linear, respectively)
! on the the method-grid (x_coords, y_coords, z_coords)
! for the subgrid coords sent by the destination process.
!
! Subroutine "PSMILe_Neigh_tricu_irreg2" is designed for (target) grids
! of type "PRISM_Irrlonlat_regvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created for linear interpolation
! 05.05.09    R. Redler    modified for cubic interpolation
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_tricu_irreg2.F90 2687 2010-10-28 15:15:52Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_tricu_irreg2.F90 2687 2010-10-28 15:15:52Z coquart $'
!
!----------------------------------------------------------------------
!
      data nca /n_corners_1d, n_corners_2d, n_corners_3d/
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
!===> Control interpolation mode
!
      if (interpolation_mode < 1 .or. interpolation_mode > ndim_3d) then
         ierrp (1) = interpolation_mode
         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_tricu_irreg2', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      endif
!
!===> Initialization
!
!  lenijk = Number of entries in all   3 directions
!
      ierror = 0
!
      nprev = npreva
      n_corners = nca (interpolation_mode)
!
      lenijk = nlocs (1) * nlocs (2)
      nbeg = nprev + 1
      nend = nprev + lenijk
!
      length (1:ndim_3d) = grid_valid_shape(2,1:ndim_3d) - &
                           grid_valid_shape(1,1:ndim_3d) + 1
!
#ifdef PRISM_ASSERTION
!
!===> Is the dimension of array "neighbors_3d" sufficient
!
      if (num_neigh < n_corners) then
         print 9970, trim(ch_id), num_neigh, n_corners
         call psmile_assert (__FILE__, __LINE__, &
                             'Number of neighbors is insufficient')
      endif
!
!===> 
!
      if (nloc < nprev + nlocs (1) * nlocs (2) ) then
         print *, trim(ch_id), " nloc, nprev, nlocs ", nloc, nprev, nlocs
         call psmile_assert (__FILE__, __LINE__, &
                             'nloc < nprev + PRODUCT (nlocs) ')
      endif
!
!===> Are the locations within the correct shape ?
!
!cdir vector
         do i = 1, nlocs(1)
            if ( srcloc(1,i) < grid_valid_shape(1,1) -1 .or. &
                 srcloc(1,i) > grid_valid_shape(2,1)    .or. &
                 srcloc(2,i) < grid_valid_shape(1,2) -1 .or. &
                 srcloc(2,i) > grid_valid_shape(2,2) ) exit
         end do

         if (i <= nlocs(1)) then
            print *, "grid valid shape is ", grid_valid_shape
            print *, "Incorrect index in srcloc, i", i, srcloc(1,i), srcloc(2,i)
            call psmile_assert (__FILE__, __LINE__, 'wrong index')
         endif
!
!===> Are the z-locations within the correct shape ?
!
!cdir vector
         do i = 1, nlocs (2)
            if ( srclocz(i) < grid_valid_shape(1,3)-1 .or. &
                 srclocz(i) > grid_valid_shape(2,3)) exit
         end do
         if (i <= nlocs(2)) then
            print *, "Incorrect index in srclocz(i), i", i, srclocz(i)
            call psmile_assert (__FILE__, __LINE__, 'wrong index')
         endif
#endif
!
!===> Control special cases (2d hyperplanes)
!
         do i = 1, ndim_3d
         if (length(i) == 1) then
            neighbors_3d (i, nbeg:nend, 1:n_corners) = grid_valid_shape(1,i)
         endif
         end do
!
!===> Generate entries (assuming mask is not available or all
!     mask points are defined)
!
         neighbors_3d (1:2, nprev+1:nprev+nlocs(1), 6) = srcloc (1:2, 1:nlocs(1))
!
!===> Duplicate first 2 indices for third direction
!
         do k = 2, nlocs (2)
            neighbors_3d (1:2, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 6) = &
            neighbors_3d (1:2, nprev+               1:nprev+  nlocs(1), 6)
         end do
!
!===> 
!
         do k = 1, nlocs (2)
            neighbors_3d (3, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 6) = srclocz (k)     
         end do
!
!===> ... I indices of neighbour indices
!
         if (length(1) > 1) then
            neighbors_3d (1, nbeg:nend, 2) = neighbors_3d (1, nbeg:nend, 6) + 1

            do n = 0, 3
               neighbors_3d (1, nbeg:nend, 1+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) - 1
               neighbors_3d (1, nbeg:nend, 2+(n*4) ) = neighbors_3d (1, nbeg:nend, 6)
               neighbors_3d (1, nbeg:nend, 3+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) + 1
               neighbors_3d (1, nbeg:nend, 4+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) + 2
            enddo

            if (n_corners > 16) &
                 neighbors_3d (1, nbeg:nend, 17:32) = neighbors_3d (1, nbeg:nend,  1:16)

         endif
!
!===> ... J indices of neighbour indices
!
         if (length(2) > 1) then

            do n = 1, 4
               neighbors_3d (2, nbeg:nend,    n) = neighbors_3d (2, nbeg:nend, 6) - 1
               neighbors_3d (2, nbeg:nend,  4+n) = neighbors_3d (2, nbeg:nend, 6)
               neighbors_3d (2, nbeg:nend,  8+n) = neighbors_3d (2, nbeg:nend, 6) + 1
               neighbors_3d (2, nbeg:nend, 12+n) = neighbors_3d (2, nbeg:nend, 6) + 2
            end do

            if (n_corners > 16) then
               neighbors_3d (2, nbeg:nend, 17:20) = neighbors_3d (2, nbeg:nend,  1: 4)
               neighbors_3d (2, nbeg:nend, 21:24) = neighbors_3d (2, nbeg:nend,  5: 8)
               neighbors_3d (2, nbeg:nend, 25:28) = neighbors_3d (2, nbeg:nend,  9:12)
               neighbors_3d (2, nbeg:nend, 29:32) = neighbors_3d (2, nbeg:nend, 13:16)
            endif

         endif
!
!===> ... K indices of neighbour indices
!
         if (length(3) > 1) then

            do n = 1, 16
               neighbors_3d (3, nbeg:nend, n) = neighbors_3d (3, nbeg:nend, 6)
            enddo

            if (n_corners > 16) then
               do n = 17, 32
                  neighbors_3d (3, nbeg:nend,n) = neighbors_3d (3, nbeg:nend, 1) + 1
               enddo
            endif

         else

            neighbors_3d (3, nbeg:nend,2:16) = srclocz (1)

         endif
!
!===> ... set indices for cyclic coordinate directions
!
         do j = 1, ndim_3d
         if (cyclic(j) .and. length(j) > 1) then
!
               do n = 1, n_corners
!cdir vector
                  do i = nbeg, nend
                     
                  if  ( neighbors_3d (j, i, n) > grid_valid_shape(2,j) ) then

                        neighbors_3d (j, i, n) = &
                        neighbors_3d (j, i, n) - length(j)

                  else if ( neighbors_3d (j, i, n) < grid_valid_shape(1,j) ) then

                        neighbors_3d (j, i, n) = &
                        neighbors_3d (j, i, n) + length(j)

                  endif
                  end do
               end do ! n
         endif
         end do ! j
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id)

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

#ifdef PRISM_ASSERTION
9970 format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
#endif

      end subroutine psmile_neigh_tricu_irreg2