psmile_neigh_tricu_3d.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_3d
!
! !INTERFACE:

      subroutine psmile_neigh_tricu_3d (                     &
               grid_valid_shape, interpolation_mode, cyclic, &
               srcloc, nloc, neighbors_3d, num_neigh, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_tricu_3d

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: nloc
!
!     Total number of locations to be transferred
!
      Integer, Intent (In)            :: srcloc (ndim_3d, nloc)
!
!     Indices of the (method) grid cell in which the point was found.

      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 64 corners)
!        interpolation_mode = 2 : Bilinear or bicubic interpolation
!                                 (determines 16 corners)
!        interpolation_mode = 1 : Linear   interpolation
!                                 (determines 4 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_3d;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  N_CORNERS_3D = Number of corners for a 3-d interpolation
!               = 4 ** ndim_3d
!  N_CORNERS_2D = Number of corners for a 2-d interpolation
!               = 4 ** ndim_2d
!
      Integer, parameter              :: n_corners_3d = 4 * 4 * 4
      Integer, parameter              :: n_corners_2d = 4 * 4
      Integer, parameter              :: n_corners_1d = 4
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: i, j, n
      Integer                         :: n_corners
!
!     ... For locations controlled
!
      Integer, Save                   :: ijkstd (ndim_3d, 1:n_corners_3d)
      Integer                         :: ijkctl (ndim_3d, 1:n_corners_3d)
      Integer                         :: length (ndim_3d)
      Integer                         :: nca (ndim_3d)
!
!     ... for error handling
!
      Integer, parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_tricu_3d" 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_3d" is designed for (target) grids
! of type "PRISM_Irrlonlatvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created for trili
! 05.05.09    R. Redler    modified for cubic interpolation
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_tricu_3d.F90 2687 2010-10-28 15:15:52Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_tricu_3d.F90 2687 2010-10-28 15:15:52Z coquart $'
!
!----------------------------------------------------------------------
!
!  Relative indices of corner (Fortran ordering)
!
      data ((ijkstd (i, n), i=1,ndim_3d), n = 1, n_corners_3d)            &
!
                            /  -1,-1, 0,   0,-1, 0,   1,-1, 0,   2,-1, 0, &
                               -1, 0, 0,   0, 0, 0,   1, 0, 0,   2, 0, 0, &
                               -1, 1, 0,   0, 1, 0,   1, 1, 0,   2, 1, 0, &
                               -1, 2, 0,   0, 2, 0,   1, 2, 0,   2, 2, 0, &
!
                               -1,-1,-1,   0,-1,-1,   1,-1,-1,   2,-1,-1, &
                               -1, 0,-1,   0, 0,-1,   1, 0,-1,   2, 0,-1, &
                               -1, 1,-1,   0, 1,-1,   1, 1,-1,   2, 1,-1, &
                               -1, 2,-1,   0, 2,-1,   1, 2,-1,   2, 2,-1, &
!
                               -1,-1, 1,   0,-1, 1,   1,-1, 1,   2,-1, 1, &
                               -1, 0, 1,   0, 0, 1,   1, 0, 1,   2, 0, 1, &
                               -1, 1, 1,   0, 1, 1,   1, 1, 1,   2, 1, 1, &
                               -1, 2, 1,   0, 2, 1,   1, 2, 1,   2, 2, 1, &
!
                               -1,-1, 2,   0,-1, 2,   1,-1, 2,   2,-1, 2, &
                               -1, 0, 2,   0, 0, 2,   1, 0, 2,   2, 0, 2, &
                               -1, 1, 2,   0, 1, 2,   1, 1, 2,   2, 1, 2, &
                               -1, 2, 2,   0, 2, 2,   1, 2, 2,   2, 2, 2 /
!
      data nca /n_corners_1d, n_corners_2d, n_corners_3d/
!
!----------------------------------------------------------------------
!
!===> Prologue
!
#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_3d', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      endif
!
!===> Initialization
!
      ierror  = 0
!
      n_corners = nca (interpolation_mode)
!
      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
!
!===> Are the locations within the correct shape ?
!     Note: the scrloc's are source locations in the 
!           method grid which has an virtual cell 0.
!
!cdir vector
         do i = 1, nloc
!
         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)   .or. &
             srcloc(3,i) < grid_valid_shape(1,3)-1 .or. &
             srcloc(3,i) > grid_valid_shape(2,3)) exit
         end do
!
      if (i <= nloc) then
         print *, "Incorrect index in i", i, srcloc (:, i)
         print *, grid_valid_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index')
      endif
#endif
!
!===> For special cases (2d Hyperplanes)
      ijkctl = ijkstd
!
         do j = 1, ndim_3d
         if (length(j) == 1) ijkctl (j, :) = 0
         end do
!
!===> Generate entries (assuming mask is not available or all
!     mask points are defined)
!
      neighbors_3d (:, 1:nloc, 6) = srcloc (:, 1:nloc)
!
         do n = 1, n_corners
         if (n == 6) cycle
!cdir vector
            do i = 1, nloc
            neighbors_3d (:, i, n) = neighbors_3d (:, i, 6) + ijkctl (:, n)
            end do
         end do
!
!===> ... 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 = 1, nloc
                  if (neighbors_3d (j, i, n) < grid_valid_shape(1,j)) then
                      neighbors_3d (j, i, n) = &
                      neighbors_3d (j, i, n) + length (j)

                  else if ( neighbors_3d (j, i, n) > grid_valid_shape(2,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_3d')
9980 format (1x, a, ': psmile_neigh_tricu_3d: eof', i6)

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

      end subroutine psmile_neigh_tricu_3d