psmile_neigh_trili_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_trili_3d
!
! !INTERFACE:

      subroutine psmile_neigh_trili_3d (                     &
               grid_valid_shape, interpolation_mode, cyclic, &
               srcloc, nloc, neighbors_3d, num_neigh, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_trili_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 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:
!
      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_trili_3d;
!             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
      Integer, Parameter              :: n_corners_2d = 2 * 2
      Integer, Parameter              :: n_corners_1d = 2
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: i, j, n
      Integer                         :: n_corners
!
!     ... For locations controlled
!
      Integer, Save                   :: ijkstd (ndim_3d, 2:n_corners_3d)
      Integer                         :: ijkctl (ndim_3d, 2: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_trili_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_trili_3d" is designed for (source) grids
! of type "PRISM_Irrlonlatvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_trili_3d.F90 2687 2010-10-28 15:15:52Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_trili_3d.F90 2687 2010-10-28 15:15:52Z coquart $'
!
!----------------------------------------------------------------------
!
!  Relative indices of corner
!
#ifdef FORTRAN_ORDER
      data ((ijkstd (i, n), i=1,ndim_3d), n = 2, n_corners_3d) &
                            /             1, 0, 0, &
                               0, 1, 0,   1, 1, 0, &
                               0, 0, 1,   1, 0, 1, &
                               0, 1, 1,   1, 1, 1/

#else
!huhu Offenbar erwartet Damien's Routine eine andere Ordnung
      data ((ijkstd (i, n), i=1,ndim_3d), n = 2, n_corners_3d) &
                            /             1, 0, 0, &
                               1, 1, 0,   0, 1, 0, &
                               0, 0, 1,   1, 0, 1, &
                               1, 1, 1,   0, 1, 1/
#endif
!
      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_trili', &
                             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)
!     ??? Braucht man hier nicht ein spezielles Flag,
!         dass das globale Gebiet eine Hyperplane ist ?
!
      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)
!
#define TRANSF_SUP_2D_INTERP
!
#ifdef TRANSF_SUP_2D_INTERP
!
! Note: The transformer doesn't currently perform 2d interpolation
!       for a point which lies on the face of the "last cell"

      neighbors_3d (:, 1:nloc, 1) = srcloc (:, 1:nloc)
#else

      if ( n_corners == N_CORNERS_2D ) then
!
!           ... 2-dimensional interpolation
!
            do j = 1, ndim_2d
            if (cyclic (j) .or. length(j) == 1) then
               neighbors_3d (j, :, 1) = srcloc (j, :)
            else
!cdir vector
               do i = 1, nloc
                  neighbors_3d (j, i, 1) = min (srcloc (j, i), &
                     grid_valid_shape(2,j)-1)
               end do

            endif
            end do ! j

!           Index in z-direction
            
         neighbors_3d (3, :, 1) = srcloc (3, :)

      else

            do j = 1, ndim_3d
            if (cyclic (j) .or. length(j) == 1) then
               neighbors_3d (j, :, 1) = srcloc (j, :)
            else
!cdir vector
               do i = 1, nloc
                  neighbors_3d (j, i, 1) = min (srcloc (j, i), &
                       grid_valid_shape(2,j)-1)
               end do
            endif
            end do

      endif

#endif /* TRANSF_SUP_2D_INTERP */
!
         do n = 2, n_corners
!cdir vector
            do i = 1, nloc
!           neighbors_3d (:, i, n) = srcloc (:, i) + ijkctl (:, n)
            neighbors_3d (:, i, n) = neighbors_3d (:, i, 1) + 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_trili_3d')
9980 format (1x, a, ': psmile_neigh_trili_3d: eof')

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

      end subroutine psmile_neigh_trili_3d