psmile_neigh_tricu_gauss2_irreg.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_irreg
!
! !INTERFACE:

      subroutine psmile_neigh_tricu_gauss2_irreg (                &
               grid_id, grid_valid_shape, interpolation_mode,     &
               srcloc, srclocz, virtual_cell, nlocs, nloc, nprev, &
               virtual_cell_available, neighbors_3d, num_neigh,   &
               neigh_bascule, ierror)
!
! !USES:
!
      use PRISM_constants
      use PSMILe, dummy_interface => PSMILe_Neigh_tricu_gauss2_irreg
      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)            :: 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)            :: nprev
!
!     Number of locations for which neighbours were already searched.
!
      Integer, Intent (In)            :: virtual_cell (nlocs(1))

!     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-cubic interpolation
!                                 (determines 64 corners) not supported
!        interpolation_mode = 2 : Bilcubic interpolation
!                                 (determines 16 corners)
!
      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_irreg;
!             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_irreg" 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"
!  and 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
! 31.01.11    M. Hanke     rewritten based on psmile_neigh_tricu_gauss2.F90
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_tricu_gauss2_irreg.F90 3032 2011-03-18 17:29:36Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_tricu_gauss2_irreg.F90 3032 2011-03-18 17:29:36Z hanke $'
!
!----------------------------------------------------------------------

#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_irreg', &
                             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
!
!===>
!
      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 ( abs(srcloc(1,i)) < grid_valid_shape(1,1) .or. &
                 abs(srcloc(1,i)) > grid_valid_shape(2,1) ) 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)
            call psmile_assert (__FILE__, __LINE__, 'wrong index')
         endif

!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 srclocs(2), i", i, srclocz(i), &
                     grid_valid_shape(1,3)-1, grid_valid_shape(2,3)
            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 i = 1, nlocs(1)

            n = i + nprev

            if (virtual_cell(i) == 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,i))))

            else

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

            endif
         enddo ! í
      else
         do i = 1, nlocs(1)

            n = i + nprev

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

#ifdef DEBUG_TRACE
      if (ictl > nprev .and. ictl <= nprev + nlocs(1)) 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,nprev+1:nprev+nlocs(1),6) >= &
             size(grids(grid_id)%reduced_gauss_data%nbr_points_per_lat)-1 .or. &
             neighbors_3d(2,nprev+1:nprev+nlocs(1),6) == 1)

         neigh_bascule(nprev+1:nprev+nlocs(1)) = 1

      else where

         neigh_bascule(nprev+1:nprev+nlocs(1)) = 0

      end where

#ifdef DEBUG_TRACE
      if (ictl > nprev .and. ictl <= nprev + nlocs(1)) 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,nprev+1:nprev+nlocs(1),i) = &
            psmile_gauss_3d_to_global_1d (grid_id, neighbors_3d(:,nprev+1:nprev+nlocs(1),i), nlocs(1))
         neighbors_3d(2,nprev+1:nprev+nlocs(1),i) = 1
#ifdef DEBUG_TRACE
         if (ictl > nprev .and. ictl <= nprev + nlocs(1)) then
            ictl_neigh_global_1d(i) = neighbors_3d (1, ictl, i)
         endif
#endif

         neighbors_3d(1,nprev+1:nprev+nlocs(1),i) = &
            psmile_gauss_1d_global_to_local (grid_id, neighbors_3d(1,nprev+1:nprev+nlocs(1),i), nlocs(1), 0)

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

#ifdef DEBUG_TRACE
      if (ictl > nprev .and. ictl <= nprev + nlocs(1)) 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 n = 1 + nprev, nprev + nlocs(1)
            neighbors_3d (1:2, n, 17:32) = neighbors_3d (1:2, n, 1:16)
         enddo
      endif
!
!===> Duplicate first two index for third direction
!
      do i = 2, nlocs (2)
         neighbors_3d (1:2, nprev+(i-1)*nlocs(1)+1:nprev+i*nlocs(1),:) = &
         neighbors_3d (1:2, nprev+               1:nprev+  nlocs(1),:)
      end do
!
!===> Set third index
!
      do i = 1, nlocs (2)
         neighbors_3d (3, nprev+(i-1)*nlocs(1)+1:nprev+i*nlocs(1), :) = srclocz (i)
      end do
!
!===> ... K indices of neighbour indices
!
      if (grid_valid_shape(2,3) - grid_valid_shape(1,3) > 0 &
          .and. num_neigh == 32) then
         do n = 17, 32
               neighbors_3d (3, nprev+1:nprev+nlocs(1)*nlocs(2),n) = &
                  neighbors_3d (3, nprev+1:nprev+nlocs(1)*nlocs(2), 1) + 1
         enddo
      endif
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id)

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

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

      end subroutine psmile_neigh_tricu_gauss2_irreg