psmile_neigh_trili_gauss2.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_Neigh_trili_gauss2
!
! !INTERFACE:

      subroutine psmile_neigh_trili_gauss2 (                       &
               grid_id, grid_valid_shape, interpolation_mode,      &
               srcloc, virtual_cell, nloc, virtual_cell_available, &
               neighbors_3d, num_neigh, ierror)
!
! !USES:
!
      use PRISM_constants
      use PSMILe, dummy_interface => PSMILe_Neigh_trili_gauss2
      use psmile_grid_reduced_gauss
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      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)            :: virtual_cell (nloc)

!     Code whether destination point is located in virtual cell with
!
!     virtual_cell (n) = 0 : Not located in a virtual cell
!     virtual_cell (n) > 0 :     Located in a virtual cell
!                      = 1 : (-1,-1) direction
!                      = 2 : ( 0,-1) direction
!                      = 3 : ( 1,-1) direction
!                      = 4 : (-1, 0) direction
!                      = 6 : ( 1, 0) direction
!                      = 7 : (-1, 1) direction
!                      = 8 : ( 0, 1) direction
!                      = 9 : ( 1, 1) direction

      Integer, Intent (In)            :: grid_id

!     Grid handle

      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) not supported
!
      Integer,            Intent (In) :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d"
!
      Logical,            Intent (In) :: virtual_cell_available
!
!     Is the virtual cell info available ?
!
! !OUTPUT PARAMETERS:
!
!     Location of indices for neighbors_3d
!
!     ^
!     |  4    3
!
!     j  1    2
!
!       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_tril_gauss2;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
      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/))
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: i, n
      Integer                         :: n_corners

#ifdef DEBUG_TRACE
!
!     ... for debugging
!
      integer                         :: ictl_neigh_global_1d(4), ictl_neigh_local_1d(4)
#endif
!
! !DESCRIPTION:
!
! Subroutine "psmile_neigh_tricu_gauss2" generates the stencil for the cubic
! interpolation and stores the local 1d coordinates the stencil in neighbors_3d.
!
! The subroutine is designed for (source) grids of type "PRISM_Gaussreduced_regvert".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
! 12.01.11    M. Hanke     revisited
! 18.02.11    M. Hanke     rewritten
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_trili_gauss2.F90 3011 2011-03-10 17:50:02Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_trili_gauss2.F90 3011 2011-03-10 17:50:02Z hanke $'
!
!----------------------------------------------------------------------
!
!===> Prologue
!
#ifdef VERBOSE
      print 9990, trim(ch_id)
      call psmile_flushstd
#endif /* VERBOSE */

!
!===> Initialization
!
      ierror  = 0
      if (interpolation_mode == 2) then
         n_corners = 4
      else if (interpolation_mode == 3) then
         n_corners = 8
      else
         ierror = PRISM_Error_Internal
         call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_trili', &
                             (/interpolation_mode/), 1, &
                             __FILE__, __LINE__ )
      endif
!
!===> inital checks
!
#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 ?
!
      do i = 1, nloc

         if (srcloc(1,i) < grid_valid_shape(1,1)   .or. &
             srcloc(1,i) > grid_valid_shape(2,1)   .or. &
             srcloc(2,i) < grid_valid_shape(1,2)   .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)) then

            print *, "Incorrect index in i", i, "src_loc(:,i)", srcloc (:, i)
            print *, "grid_valid_shape", grid_valid_shape
            call psmile_assert (__FILE__, __LINE__, &
                                'wrong index')
         endif
      end do
#endif

!
!==> Actual work
!
      ! 1st generate the 3d stencil

      neighbors_3d(:,:,1) = psmile_gauss_local_1d_to_3d (grid_id, abs(srcloc(1,:)), nloc)

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

            if (virtual_cell(n) == 0) then

               neighbors_3d(:,n,1:4) = psmile_gauss_get_bili_stencil (grid_id, neighbors_3d(:,n,1))

            else

               neighbors_3d(:,n,1:4) = psmile_gauss_shift_bili_stencil (grid_id, neighbors_3d(:,n,1), &
                                                                        shift_3d(:,virtual_cell(n)))

            endif
         enddo ! n
      else
         do n = 1, nloc
            neighbors_3d(:,n,1:4) = psmile_gauss_get_bili_stencil (grid_id, neighbors_3d(:,n,1))
         enddo ! n
      endif

#ifdef DEBUG_TRACE
      if (ictl > 0 .and. ictl <= nloc) then
         print *, "neighbors_3d (3d):"
         do i = 1, 4
            print *, ' neighbors_3d (:,ictl,i) ', i, ':', neighbors_3d (:, ictl, i)
         enddo
      endif
#endif

      ! 2nd transform the 3d coordinates into local 1d ones

      do i = 1, 4

         neighbors_3d(1,:,i) = psmile_gauss_3d_to_global_1d (grid_id, neighbors_3d(:,:,i), nloc)
         neighbors_3d(2:3,:,i) = 1
#ifdef DEBUG_TRACE
         if (ictl > 0 .and. ictl < nloc) then
            ictl_neigh_global_1d(i) = neighbors_3d (1, ictl, i)
         endif
#endif

         neighbors_3d(1,:,i) = psmile_gauss_1d_global_to_local (grid_id, neighbors_3d(1,:,i), nloc, 0)

#ifdef DEBUG_TRACE
         if (ictl > 0 .and. ictl < nloc) then
            ictl_neigh_local_1d(i) = neighbors_3d (1, ictl, i)
         endif
#endif
      enddo

#ifdef DEBUG_TRACE
      if (ictl > 0 .and. ictl < nloc) then
         print *, "neighbors_3d (1d global):"
         print *, ' neighbors_3d (1,ictl,1:4) ', ictl_neigh_global_1d
         print *, "neighbors_3d (1d local):"
         print *, ' neighbors_3d (1,ictl,1:4) ', ictl_neigh_local_1d
      endif
#endif

!
      if ( interpolation_mode == ndim_3d ) then
!cdir vector
         do i = 1, nloc
            neighbors_3d (1:2, i, 5:8) = neighbors_3d (1:2, i, 1:4)
         enddo
      endif
!
!===> ... K indices of neighbour indices
!
      do i = 1, n_corners
         neighbors_3d (3, 1:nloc, i) = srcloc(3,1:nloc)
      enddo

      if (interpolation_mode == ndim_3d .and. &
          grid_valid_shape(2,3) - grid_valid_shape(1,3) > 0) then

         do i = 5, 8
            neighbors_3d (3, 1:nloc, i) = neighbors_3d (3, 1:nloc, i) + 1
         enddo
      endif
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id)

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

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

      end subroutine psmile_neigh_trili_gauss2