psmile_neigh_tricu_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_tricu_gauss2
!
! !INTERFACE:

      subroutine psmile_neigh_tricu_gauss2 (                       &
               grid_id, grid_valid_shape, interpolation_mode,      &
               srcloc, virtual_cell, nloc, virtual_cell_available, &
               neighbors_3d, num_neigh, neigh_bascule, ierror)
!
! !USES:
!
      use PRISM_constants
      use PSMILe, dummy_interface => PSMILe_Neigh_tricu_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 64 corners) not supported
!        interpolation_mode = 2 : Bilinear or bicubic interpolation
!                                 (determines 16 corners)
!        interpolation_mode = 1 : Linear   interpolation
!                                 (determines 4 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
!
!           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)           :: neigh_bascule(nloc)

!     Indicator for special index ordering

      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_Neigh_tricu_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
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)

#ifdef DEBUG_TRACE
!
!     ... for debugging
!
      integer                         :: ictl_neigh_global_1d(16), ictl_neigh_local_1d(16)
#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.
! In addition it stores the information whether a stencil covers the pole or
! not in neigh_bascule.
!
! The subroutine is designed for (source) grids of type "PRISM_Gaussreduced_regvert".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created for trili
! 05.05.09    R. Redler    modified for cubic interpolation
! 27.01.11    M. Hanke     rewritten
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_tricu_gauss2.F90 2968 2011-02-18 15:26:34Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_tricu_gauss2.F90 2968 2011-02-18 15:26:34Z hanke $'
!
!----------------------------------------------------------------------
!
!===> Prologue
!
#ifdef VERBOSE
      print 9990, trim(ch_id)
      call psmile_flushstd
#endif /* VERBOSE */


!
!===> inital checks
!
#ifdef PRISM_ASSERTION
!
!===> Control interpolation mode
!
      if (interpolation_mode /= 2) then
         ierrp (1) = interpolation_mode
         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_tricu_gauss2', &
                             ierrp, 1, __FILE__, __LINE__ )
         return
      endif
!
!===> Is the dimension of array "neighbors_3d" sufficient
!
      if (num_neigh < 16) then
         print 9970, trim(ch_id), num_neigh, 16
         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
!           (negative for Gauss Grid).
!
      do i = 1, nloc
!
         if (abs(srcloc(1,i)) < grid_valid_shape(1,1)   .or. &
             abs(srcloc(1,i)) > grid_valid_shape(2,1)   .or. &
             abs(srcloc(2,i)) < grid_valid_shape(1,2)   .or. &
             abs(srcloc(2,i)) > grid_valid_shape(2,2)   .or. &
             abs(srcloc(3,i)) < grid_valid_shape(1,3)-1 .or. &
             abs(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
!
!===> Initialization
!
      ierror  = 0

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

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

            if (virtual_cell(n) == 0) then

               neighbors_3d(:,n,1:16) = psmile_gauss_get_bicu_stencil (grid_id, &
                                          psmile_gauss_local_1d_to_3d (grid_id, abs(srcloc(1,n))))

            else

               neighbors_3d(:,n,1:16) = psmile_gauss_shift_bicu_stencil (grid_id,                  &
                                          psmile_gauss_local_1d_to_3d (grid_id, abs(srcloc(1,n))), &
                                          shift_3d(:,virtual_cell(n)))

            endif
         enddo ! n
      else
         do n = 1, nloc
            neighbors_3d(:,n,1:16) = psmile_gauss_get_bicu_stencil (grid_id, &
                                       psmile_gauss_local_1d_to_3d (grid_id, abs(srcloc(1,n))))
         enddo ! n
      endif

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

      ! 2nd check for stencils that go over the pole

      where (neighbors_3d(2,:,6) >= grids(grid_id)%nbr_latitudes-1 .or. &
             neighbors_3d(2,:,6) == 1)

         neigh_bascule(:) = 1

      else where

         neigh_bascule(:) = 0

      end where

#ifdef DEBUG_TRACE
      if (ictl > 0 .and. ictl < nloc) then
         print *, "neigh_bascule(ictl)", neigh_bascule(ictl)
      endif
#endif

      ! 3rd transform the 3d coordinates into local 1d ones

      do i = 1, 16

         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:16) ', ictl_neigh_global_1d
         print *, "neighbors_3d (1d local):"
         print *, ' neighbors_3d (1,ictl,1:16) ', ictl_neigh_local_1d
      endif
#endif

!
      if ( num_neigh == 32 ) then
!cdir vector
         do i = 1, nloc
            neighbors_3d (1:2, i, 17:32) = neighbors_3d (1:2, i, 1:16)
         enddo
      endif
!
!===> ... K indices of neighbour indices
!
      do i = 1, num_neigh
         neighbors_3d (3, 1:nloc, i) = srcloc(3,1:nloc)
      enddo

      if (grid_valid_shape(2,3) - grid_valid_shape(1,3) > 0) then

         if (num_neigh == 32) then
            do i = 17, 32
               neighbors_3d (3, 1:nloc, i) = neighbors_3d (3, 1:nloc, i) + 1
            enddo
         endif
      endif
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id)

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

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

      end subroutine psmile_neigh_tricu_gauss2