psmile_reduced_gauss_utils.F90

Go to the documentation of this file.

!------------------------------------------------------------------------
! Copyright 2011, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!------------------------------------------------------------------------
!
! !DESCRIPTION:
!
! This file contains utility routines and functions for reduced gauss
! grids.
! psmile_init_gauss_data: initilises the data for reduced gauss grids,
!                         that is stored in grids(see psmile.F90)
! psmile_gauss_gen_aux_grid : generates a reglonlatvrt auxiliary grid
!                             which is used in the multigrid search
!                             instead of the actual reduced gauss grid
! psmile_gauss_gen_aux_grid_map : generates a mapping between the
!                                 auxiliary grid and the reduced gauss
!                                 grid. this is used to translate the
!                                 results from the multigrid search into
!                                 reduced gauss grid coordinates
! psmile_free_aux_grid_data : frees the memory consumed by the auxiliary
!                             grid and the associated map
! psmile_gauss_blockidx_to_glat: returns the global latitude of the
!                                provided block index
! psmile_gauss_local_*_1d_to_3d: returns a global 3d index that corresponds
!                                to the provided local 1d index
! get_gauss_neighbour_cell: returns the longitude index of cell that is
!                           a neighbour to a given second longitude index.
!                           (This is used in cases where the longitude rows
!                           of both points have different numbers of
!                           points)
! get_gauss_opposite_neighbour_cell: same as get_gauss_neighbour_cell
!                                    except that this one is used in
!                                    cases where both points are on
!                                    different sides of a pole
! psmile_gauss_get_bicu_stencil: generates a interpolation stencil for bicubic
!                                interpolation based on a global 3d start
!                                location
! psmile_gauss_shift_bicu_stencil: nearly the same as
!                                  psmile_gauss_get_bicu_stencil except,
!                                  that the stencil return is shifted
!                                  according to the provided shift-values
! psmile_gauss_get_neigh_stencil: generates a stencil which contains the
!                                 point provided by the user and its 8
!                                 direct neighbours
! psmile_gauss_get_bili_stencil: generates a interpolation stencil for bilinear
!                                interpolation based on a global 3d start
!                                location
! psmile_gauss_shift_bili_stencil: nearly the same as
!                                  psmile_gauss_get_bili_stencil except,
!                                  that the stencil return is shifted
!                                  according to the provided shift-values
! psmile_gauss_3d_to_global_1d: transforms an array of reduced gauss 3d
!                               coordinates into global 1d ones
! psmile_gauss_1d_global_to_local: transforms reduced gauss global 1d
!                                  coordinates into local ones. it also
!                                  takes remote points that stored
!                                  locally ("remote_index"-array) into
!                                  account. if a global value has no
!                                  representation in the local 1d index
!                                  space it is set to the fill value
!                                  provided by the user.
! psmile_gauss_3d_to_local_1d: transforms reduced gauss 3d coordinates into
!                              local 1d ones. it also takes remote points that
!                              stored locally ("remote_index"-array) into
!                              account. if a 3d value has no representation
!                              in the local 1d index space it is set to the fill
!                              value provided by the user.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 25.01.11    M. Hanke     created
!
!----------------------------------------------------------------------
!
!  $Id: psmile_reduced_gauss_utils.F90 3021 2011-03-17 09:26:03Z hanke $
!  $Author: hanke $
!
!----------------------------------------------------------------------

subroutine psmile_init_gauss_data (grid_id, nbr_points_per_lat)

   use psmile, only : grids

   implicit none

   integer, intent(in) :: grid_id
   integer, intent(in) :: nbr_points_per_lat(:)

   integer :: i

   allocate (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat(size (nbr_points_per_lat)), &
             grids(grid_id)%reduced_gauss_data%global_lat_offsets(size (nbr_points_per_lat)))

   grids(grid_id)%reduced_gauss_data%nbr_points_per_lat = nbr_points_per_lat

   grids(grid_id)%reduced_gauss_data%global_lat_offsets(1) = 0

   do i = 2, size (nbr_points_per_lat)
      grids(grid_id)%reduced_gauss_data%global_lat_offsets(i) = &
         grids(grid_id)%reduced_gauss_data%global_lat_offsets(i-1) + nbr_points_per_lat(i-1)
   enddo

end subroutine psmile_init_gauss_data

!--------------------------------------------------------------------------------------------------

subroutine psmile_gauss_gen_aux_grid (grid_id)

   use psmile_common, only : ndim_3d
   use psmile, only : grids, MPI_REAL, psmile_assert

   implicit none

   integer, intent(in) :: grid_id

   double precision :: max_size_grid_cell (ndim_3d)
   double precision :: local_partition_extent (2, ndim_3d)

   integer :: aux_grid_size (ndim_3d)

   double precision :: dxyz(ndim_3d)

   double precision, pointer :: k_corners_dble(:)
   real, pointer             :: k_corners_real(:)

   integer :: i, j

   ! determine the resolution required for the auxiliary grid
   max_size_grid_cell(1) = 360.d0 / &
      dble (maxval (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat))

   max_size_grid_cell(2) = 180.d0 / &
      dble ( size (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat))

   if (grids (grid_id)%corner_pointer%corner_datatype == MPI_REAL) then

      k_corners_real => grids (grid_id)%corner_pointer%corners_real(3)%vector
      max_size_grid_cell(3) = minval ( abs (k_corners_real(size(k_corners_real)/2+1:size(k_corners_real)) - &
                                            k_corners_real(                       1:size(k_corners_real)/2)))
   else

      k_corners_dble => grids (grid_id)%corner_pointer%corners_dble(3)%vector
      max_size_grid_cell(3) = minval ( abs (k_corners_dble(size(k_corners_dble)/2+1:size(k_corners_dble)) - &
                                            k_corners_dble(                       1:size(k_corners_dble)/2)))
   endif

   do i = 1, ndim_3d
      if (grids (grid_id)%corner_pointer%corner_datatype == MPI_REAL) then

         ! compute the extent of the local partition
         local_partition_extent(1, i) = minval (grids (grid_id)%corner_pointer%corners_real(i)%vector)
         local_partition_extent(2, i) = maxval (grids (grid_id)%corner_pointer%corners_real(i)%vector)
      else

         ! compute the extent of the local partition
         local_partition_extent(1, i) = minval (grids (grid_id)%corner_pointer%corners_dble(i)%vector)
         local_partition_extent(2, i) = maxval (grids (grid_id)%corner_pointer%corners_dble(i)%vector)
      endif
   enddo

#ifdef PRISM_ASSERTION
   if (any (max_size_grid_cell == 0) .or. &
       any (local_partition_extent(2,:) - &
            local_partition_extent(1,:) <= 0)) then

         print *, "grid extent in at least one dimension is equal to 0"
         call psmile_assert (__FILE__, __LINE__, &
                             'error in grid corner data')
   endif
#endif

   aux_grid_size (:) = ceiling ((local_partition_extent(2,:) - local_partition_extent(1,:)) / &
                                max_size_grid_cell(:))

   ! at this point the resolution of the auxiliary grid is increaded by a factor of 3
   ! (due to the properties of reduced gauss grids a factor of 2 would be too less)
   ! in order to make mapping of auxiliary grid cells to the actual cell of the reduced
   ! gauss grid more accurate
   where (aux_grid_size(:) > 1) ! it makes no sense to increase resoltion if there is
                                ! only one cell in the respective dimension
      aux_grid_size(:) = aux_grid_size(:) * 3
   endwhere

   ! compute the actual auxiliary grid
   allocate (grids(grid_id)%reduced_gauss_data%aux_corners(1)%vector(2*aux_grid_size(1)), &
             grids(grid_id)%reduced_gauss_data%aux_corners(2)%vector(2*aux_grid_size(2)), &
             grids(grid_id)%reduced_gauss_data%aux_corners(3)%vector(2*aux_grid_size(3)))

   dxyz(:) = (local_partition_extent(2,:) - local_partition_extent(1,:)) / aux_grid_size (:)

   do i = 1, ndim_3d
      grids(grid_id)%reduced_gauss_data%aux_corners(i)%vector(2*aux_grid_size(i)) = &
                                                                  local_partition_extent(2,i)
   enddo

   do i = 1, ndim_3d
      do j = 1, aux_grid_size(i)

         grids(grid_id)%reduced_gauss_data%aux_corners(i)%vector(j) = &
               local_partition_extent(1,i) + (j - 1) * dxyz(i)
      enddo

      grids(grid_id)%reduced_gauss_data%aux_corners(i)%vector(aux_grid_size(i)+1:2*aux_grid_size(i)-1) = &
         grids(grid_id)%reduced_gauss_data%aux_corners(i)%vector(2:aux_grid_size(i))
   enddo

end subroutine psmile_gauss_gen_aux_grid

!--------------------------------------------------------------------------------------------------

subroutine psmile_gauss_gen_aux_grid_map (grid_id)

   use psmile_common, only : ndim_3d, dble_vector, real_vector
   use psmile, only        : grids, MPI_REAL

   implicit none

   integer, intent (in) :: grid_id

   integer :: cell_index, num_total_cells, grid_extent(2)
   integer :: iblock, i, j, level
   integer :: cell_overlap_range(2,ndim_3d)

   type (dble_vector), pointer :: p_aux_corners(:), p_corners_dble(:)
   type (real_vector), pointer :: p_corners_real(:)

   double precision :: cell (2, ndim_3d)

   p_aux_corners => grids(grid_id)%reduced_gauss_data%aux_corners

   if (grids (grid_id)%corner_pointer%corner_datatype == MPI_REAL) then
      p_corners_real => grids(grid_id)%corner_pointer%corners_real
   else
      p_corners_dble => grids(grid_id)%corner_pointer%corners_dble
   endif

   allocate (grids(grid_id)%reduced_gauss_data%aux_3d_to_local_1d_map( &
               size (p_aux_corners(1)%vector)/2,                       &
               size (p_aux_corners(2)%vector)/2,                       &
               size (p_aux_corners(3)%vector)/2))

   ! initialise with "outside"
   grids(grid_id)%reduced_gauss_data%aux_3d_to_local_1d_map = -1

   cell_index = 1
   num_total_cells = sum (grids(grid_id)%extent(:,1))
   grid_extent(:) = grids(grid_id)%grid_shape(2,1:2) - grids(grid_id)%grid_shape(1,1:2) + 1

   ! grid_extent(1) should be the same as grids(grid_id)%extent

   ! for all blocks in the local partition
   do iblock = 1, size (grids(grid_id)%extent, 1)

      ! for all horizontal level
      do level = 1, grid_extent(2)

         if (grids (grid_id)%corner_pointer%corner_datatype == MPI_REAL) then
            ! get extent of cell in 3rd dimension
            cell(:, 3) = p_corners_real(3)%vector(level:level+grid_extent(2):grid_extent(2))
         else
            ! get extent of cell in 3rd dimension
            cell(:, 3) = p_corners_dble(3)%vector(level:level+grid_extent(2):grid_extent(2))
         endif

         ! determine the index range of the cell in the auxiliary grid
         ! for the 3rd dimension
         cell_overlap_range (:,3) = get_matching_range (cell(:, 3), p_aux_corners(3))

         ! for all cells in the current block (one block only contains cells with a common latitude)
         do i = 1, grids(grid_id)%extent(iblock, 1)

            do j = 1, 2

               if (grids (grid_id)%corner_pointer%corner_datatype == MPI_REAL) then
                  ! get the extent of the cell in the 1st and 2nd dimension
                  cell(:, j) = p_corners_real(j)%vector(cell_index:cell_index+num_total_cells:num_total_cells)
               else
                  ! get the extent of the cell in the 1st and 2nd dimension
                  cell(:, j) = p_corners_dble(j)%vector(cell_index:cell_index+num_total_cells:num_total_cells)
               endif

               ! determine the index range of the cell in the auxiliary grid
               ! for the 1st and 2nd dimension
               cell_overlap_range (:,j) = get_matching_range (cell(:, j), p_aux_corners(j))
            enddo ! j

            ! set the mapping of the auxiliary to the cell
            grids(grid_id)%reduced_gauss_data%aux_3d_to_local_1d_map (&
                     cell_overlap_range(1,1):cell_overlap_range(2,1), &
                     cell_overlap_range(1,2):cell_overlap_range(2,2), &
                     cell_overlap_range(1,3):cell_overlap_range(2,3)) = cell_index

            ! set cell_index to the next cell
            cell_index = cell_index + 1
         enddo ! i
      enddo ! level
   enddo ! iblock

contains

   function get_matching_range (cell_corners, grid_corners)

      double precision, intent (in)   :: cell_corners(2)
      type (dble_vector), intent (in) :: grid_corners

      integer :: get_matching_range(2)

      integer :: i, num_cells_in_grid

      num_cells_in_grid = size (grid_corners%vector) / 2

      ! determine the lower bound of the index range of the cell in the grid
      do i = 1, num_cells_in_grid

         if (grid_corners%vector(i) > cell_corners(1)) exit
      enddo ! i
      get_matching_range(1) = max (1, i-1)

      ! determine the upper bound of the index range of the cell in the auxiliary grid
      ! for the 3rd dimension
      do i = 1, num_cells_in_grid

         if (grid_corners%vector(i + num_cells_in_grid) >= cell_corners(2)) exit
      enddo ! i
      get_matching_range(2) = min (i, num_cells_in_grid)

   end function get_matching_range

end subroutine psmile_gauss_gen_aux_grid_map

!--------------------------------------------------------------------------------------------------

subroutine psmile_free_aux_grid_data (grid_id)

   use psmile, only : grids

   implicit none

   integer, intent (in) :: grid_id

   deallocate (grids(grid_id)%reduced_gauss_data%aux_3d_to_local_1d_map, &
               grids(grid_id)%reduced_gauss_data%aux_corners(1)%vector,  &
               grids(grid_id)%reduced_gauss_data%aux_corners(2)%vector,  &
               grids(grid_id)%reduced_gauss_data%aux_corners(3)%vector)

end subroutine psmile_free_aux_grid_data

!--------------------------------------------------------------------------------------------------

integer function psmile_gauss_blockidx_to_glat (grid_id, block_idx)

   use psmile, only : grids

   implicit none

   integer, intent (in) :: grid_id, block_idx

   integer :: block_start_offset, accu, i

   block_start_offset = grids(grid_id)%partition(block_idx, 1)

   accu = 0

   do i = 1, size (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat)

      accu = accu + grids(grid_id)%reduced_gauss_data%nbr_points_per_lat(i)

      if (accu > block_start_offset) exit
   enddo

   psmile_gauss_blockidx_to_glat = i

end function psmile_gauss_blockidx_to_glat

!--------------------------------------------------------------------------------------------------

function psmile_gauss_1_local_1d_to_3d (grid_id, idx)

   use psmile_common, only             : ndim_3d
   use psmile, only                    : grids
   use psmile_grid_reduced_gauss, only : psmile_gauss_blockidx_to_glat

   implicit none

   integer, intent (in) :: grid_id, idx

   integer              :: psmile_gauss_1_local_1d_to_3d (ndim_3d)

   integer :: lon_idx, i

   ! if there is additional information on the local blocks available we can
   ! optimise the transformation
   if (associated (grids(grid_id)%reduced_gauss_data%local_block_info)) then

      do i = 1, size (grids(grid_id)%reduced_gauss_data%local_block_info)

         if (idx >= grids(grid_id)%reduced_gauss_data%local_block_info(i)%local_1d_start_idx .and. &
             idx <= grids(grid_id)%reduced_gauss_data%local_block_info(i)%local_1d_end_idx) then

            psmile_gauss_1_local_1d_to_3d(:) = &
               grids(grid_id)%reduced_gauss_data%local_block_info(i)%global_3d_start_idx

            psmile_gauss_1_local_1d_to_3d(1) = psmile_gauss_1_local_1d_to_3d(1) + idx - &
               grids(grid_id)%reduced_gauss_data%local_block_info(i)%local_1d_start_idx

            exit
         endif
      enddo ! i

   else ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

      lon_idx = idx

      ! for all extents
      do i = 1, size (grids(grid_id)%extent, 1)

         if (lon_idx <= grids(grid_id)%extent(i,1)) exit

         lon_idx = lon_idx - grids(grid_id)%extent(i,1)
      enddo

      psmile_gauss_1_local_1d_to_3d(2) = psmile_gauss_blockidx_to_glat(grid_id, i)
      psmile_gauss_1_local_1d_to_3d(1) = lon_idx + grids(grid_id)%partition(i, 1) - &
         grids(grid_id)%reduced_gauss_data%global_lat_offsets(psmile_gauss_1_local_1d_to_3d(2))
      psmile_gauss_1_local_1d_to_3d(3) = 1

   endif ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

end function psmile_gauss_1_local_1d_to_3d

!--------------------------------------------------------------------------------------------------

function psmile_gauss_n_local_1d_to_3d (grid_id, idx, num_points)

   use psmile_common, only             : ndim_3d
   use psmile, only                    : grids
   use psmile_grid_reduced_gauss, only : psmile_gauss_blockidx_to_glat, block_info

   implicit none

   integer, intent (in)       :: num_points
   integer, intent (in)       :: grid_id, idx(num_points)
      
   integer                    :: psmile_gauss_n_local_1d_to_3d (ndim_3d,num_points)
      
   integer                    :: lon_idx, i, n, curr_block, num_blocks
   type (block_info), pointer :: local_block_info(:)


   ! if there is additional information on the local blocks available we can
   ! optimise the transformation
   if (associated (grids(grid_id)%reduced_gauss_data%local_block_info)) then

      curr_block = 1
      num_blocks = size (grids(grid_id)%reduced_gauss_data%local_block_info)
      local_block_info => grids(grid_id)%reduced_gauss_data%local_block_info

outer: do n = 1, num_points

         if (idx(n) >= local_block_info(curr_block)%local_1d_start_idx .and. &
             idx(n) <= local_block_info(curr_block)%local_1d_end_idx) then

            psmile_gauss_n_local_1d_to_3d(2:3,n) = &
               local_block_info(curr_block)%global_3d_start_idx(2:3)

            psmile_gauss_n_local_1d_to_3d(1,n) = &
               local_block_info(curr_block)%global_3d_start_idx(1) + idx(n) - &
               local_block_info(curr_block)%local_1d_start_idx

            cycle outer
         endif

         do i = curr_block+1, num_blocks, 1

            if (idx(n) >= local_block_info(i)%local_1d_start_idx .and. &
                idx(n) <= local_block_info(i)%local_1d_end_idx) then

               psmile_gauss_n_local_1d_to_3d(2:3,n) = &
                  local_block_info(i)%global_3d_start_idx(2:3)

               psmile_gauss_n_local_1d_to_3d(1,n) = &
                  local_block_info(i)%global_3d_start_idx(1) + idx(n) - &
                  local_block_info(i)%local_1d_start_idx

               curr_block = curr_block + 1
               cycle outer
            endif
         enddo ! i

         do i = curr_block-1, 1, -1

            if (idx(n) >= local_block_info(i)%local_1d_start_idx .and. &
                idx(n) <= local_block_info(i)%local_1d_end_idx) then

               psmile_gauss_n_local_1d_to_3d(2:3,n) = &
                  local_block_info(i)%global_3d_start_idx(2:3)

               psmile_gauss_n_local_1d_to_3d(1,n) = &
                  local_block_info(i)%global_3d_start_idx(1) + idx(n) - &
                  local_block_info(i)%local_1d_start_idx

               curr_block = curr_block - 1
               cycle outer
            endif
         enddo ! i
      enddo outer ! n

   else ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

      do n = 1, num_points

         lon_idx = idx(n)

         ! for all extents
         do i = 1, size (grids(grid_id)%extent, 1)

            if (lon_idx <= grids(grid_id)%extent(i,1)) exit

            lon_idx = lon_idx - grids(grid_id)%extent(i,1)
         enddo

         psmile_gauss_n_local_1d_to_3d(2,n) = psmile_gauss_blockidx_to_glat(grid_id, i)
         psmile_gauss_n_local_1d_to_3d(1,n) = lon_idx + grids(grid_id)%partition(i, 1) - &
            grids(grid_id)%reduced_gauss_data%global_lat_offsets(psmile_gauss_n_local_1d_to_3d(2,n))
         psmile_gauss_n_local_1d_to_3d(3,n) = 1
      enddo ! n

   endif ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

end function psmile_gauss_n_local_1d_to_3d

!--------------------------------------------------------------------------------------------------

integer function get_gauss_neighbour_cell (idx, nbr_points_lat, nbr_points_lat_neigh)

   implicit none

   integer, intent (in) :: idx,                ! 1-based index of the cell in the latitude band
                           nbr_points_lat,     ! number of cells in the latitude band
                           nbr_points_lat_neigh ! number of cells in the latitude band for which
                                                ! we search the neighbour cell index

   integer :: intermediate_idx

   if (nbr_points_lat == nbr_points_lat_neigh) then
      get_gauss_neighbour_cell = idx
      return
   endif

   ! Algorithm description:
   ! at first we map the points of both latitude bands to a band with
   ! (nbr_points_lat * nbr_points_lat_neigh * 3)
   ! In our intermediate latitude band we have (nbr_points_lat_neigh * 3) points for each point
   ! in nbr_points_lat. Our intermediate idx is a point in the middle of the point range that
   ! matches out given idx.
   ! In the end we compute the point in the latitude band which matches our intermediate index.

   intermediate_idx = (idx-1)*3*nbr_points_lat_neigh + nint(real(nbr_points_lat_neigh)*1.5)

   get_gauss_neighbour_cell = intermediate_idx / (3 * nbr_points_lat) + 1

end function get_gauss_neighbour_cell

!--------------------------------------------------------------------------------------------------

integer function get_gauss_opposite_neighbour_cell (idx, nbr_points_lat, nbr_points_lat_neigh)

   implicit none

   integer, intent (in) :: idx,                ! 1-based index of the cell in the latitude band
                           nbr_points_lat,     ! number of cells in the latitude band
                           nbr_points_lat_neigh ! number of cells in the latitude band for which
                                                ! we search the neighbour cell index

   integer :: intermediate_idx

   ! Algorithm description:
   ! at first we map the points of both latitude bands to a band with
   ! (nbr_points_lat * nbr_points_lat_neigh * 3)
   ! In our intermediate latitude band we have (nbr_points_lat_neigh * 3) points for each point
   ! in nbr_points_lat. Our intermediate idx is a point in the middle of the point range that
   ! matches out given idx.
   ! In the end we compute the point in the latitude band which matches our intermediate index.

   intermediate_idx = (idx-1)*3*nbr_points_lat_neigh + nint(real(nbr_points_lat_neigh)*1.5)

   ! get the intermedate index on the other side of the latitude band (turn 180 degree)
   intermediate_idx = mod(intermediate_idx + nbr_points_lat * nbr_points_lat_neigh * 3 / 2, &
                           nbr_points_lat * nbr_points_lat_neigh * 3)

   get_gauss_opposite_neighbour_cell = intermediate_idx / (3 * nbr_points_lat) + 1

end function get_gauss_opposite_neighbour_cell

!--------------------------------------------------------------------------------------------------

function psmile_gauss_get_bicu_stencil (grid_id, base)

   use psmile_common, only             : ndim_3d
   use psmile_grid_reduced_gauss, only : psmile_gauss_shift_bicu_stencil

   implicit none

   integer, intent (in) :: grid_id, base(ndim_3d)

   integer :: psmile_gauss_get_bicu_stencil(ndim_3d, 16)

   psmile_gauss_get_bicu_stencil = psmile_gauss_shift_bicu_stencil(grid_id, base, (/0,0,0/))

end function psmile_gauss_get_bicu_stencil

!--------------------------------------------------------------------------------------------------

function psmile_gauss_shift_bicu_stencil (grid_id, base, shift)

   use psmile_common, only             : ndim_3d
   use psmile, only                    : grids, psmile_assert
   use psmile_grid_reduced_gauss, only : get_gauss_opposite_neighbour_cell, &
                                         get_gauss_neighbour_cell

   implicit none

   integer, intent (in) :: grid_id, base (ndim_3d), shift(ndim_3d)

   integer :: psmile_gauss_shift_bicu_stencil (ndim_3d, 16)

   integer            :: global_num_lats, i, i_
   integer, pointer   :: g_nbr_points_per_lat(:)
   integer, parameter :: lon_shift(3) = (/1,2,-1/)
   integer            :: orientation(4)

   ! the points returned have the following relation to each other ("base" has index 6):
   !
   !    13  14  15  16
   !
   !     9  10  11  12
   !
   !     5   6   7   8
   !
   !     1   2   3   4

   g_nbr_points_per_lat => grids(grid_id)%reduced_gauss_data%nbr_points_per_lat
   global_num_lats = size (g_nbr_points_per_lat)
   orientation(:) = 1

#ifdef PRISM_ASSERTION
   if (any (abs (shift(:)) > 1)) then
      call psmile_assert (__FILE__, __LINE__, &
                          "shifting by more than one cell is not supported")
   endif
   if (base(2)+shift(2) < 1 .or. &
       base(2)+shift(2) > size (g_nbr_points_per_lat)) then
      call psmile_assert (__FILE__, __LINE__, &
                          "shifting of the base point over the pole is not supported")
   endif
#endif

   ! compute the second point

   ! if we need to recalculate the i index of the point 2
   if (shift(2) /= 1) then
   
      ! some special care needs to be taken in case it is on the other side of the pole
      if (base(2)+shift(2) == 1) then

            psmile_gauss_shift_bicu_stencil(1,2) =                                   &
               get_gauss_opposite_neighbour_cell(base(1), g_nbr_points_per_lat(1),   &
                                                          g_nbr_points_per_lat(1)) + &
               shift(1)

            psmile_gauss_shift_bicu_stencil(2,2) = 1

            orientation(1) = -1
      else

            psmile_gauss_shift_bicu_stencil(1,2) =                                &
               get_gauss_neighbour_cell(base(1),                                  &
                                        g_nbr_points_per_lat(base(2)),            &
                                        g_nbr_points_per_lat(base(2)-1+shift(2))) &
               + shift(1)

            psmile_gauss_shift_bicu_stencil(2,2) = base(2) - 1 + shift(2)
      endif
   else ! shift(2) == 1

      psmile_gauss_shift_bicu_stencil(1,2) = base(1) + shift(1)
      psmile_gauss_shift_bicu_stencil(2,2) = base(2)
   endif

   psmile_gauss_shift_bicu_stencil(3, 2) = base(3) + shift(3)

   ! compute the sixth point

   ! if we need to recalculate the i index of the point 6
   if (abs(shift(2)) == 1) then
   
      psmile_gauss_shift_bicu_stencil(1,6) =                                       &
         get_gauss_neighbour_cell(base(1), g_nbr_points_per_lat(base(2)),          &
                                           g_nbr_points_per_lat(base(2)+shift(2))) &
         + shift(1)

      psmile_gauss_shift_bicu_stencil(2,6) = base(2) + shift(2)

   else ! shift(2) == 0

      psmile_gauss_shift_bicu_stencil(1:2,6) = base(1:2) + shift(1:2)
   endif

   psmile_gauss_shift_bicu_stencil(3,6) = base(3) + shift(3)

   ! compute the eleventh point

   ! if we need to recalculate the i index of point 4
   if (shift(2) >= 0) then

      ! some special care needs to be taken in case it is on the other side of the pole
      if (base(2)+shift(2) == global_num_lats) then

            psmile_gauss_shift_bicu_stencil(1,10) =               &
               get_gauss_opposite_neighbour_cell(                 &
                  base(1), g_nbr_points_per_lat(global_num_lats), &
                           g_nbr_points_per_lat(global_num_lats)) &
               + shift(1)

            psmile_gauss_shift_bicu_stencil(2,10) = global_num_lats

            orientation(3) = -1
      else

            psmile_gauss_shift_bicu_stencil(1,10) =                               &
               get_gauss_neighbour_cell(base(1),                                  &
                                        g_nbr_points_per_lat(base(2)),            &
                                        g_nbr_points_per_lat(base(2)+1+shift(2))) &
               + shift(1)

            psmile_gauss_shift_bicu_stencil(2,10) = base(2) + 1 + shift(2)
      endif

   else ! shift(2) == -1
   
      psmile_gauss_shift_bicu_stencil(1,10) = base(1) + shift(1)
      psmile_gauss_shift_bicu_stencil(2,10) = base(2)
   endif

   psmile_gauss_shift_bicu_stencil(3,10) = base(3) + shift(3)

   ! compute the fourteenth point

   ! some special care needs to be taken in case it is on the other side of the pole
   if (base(2)+shift(2) == global_num_lats) then

         psmile_gauss_shift_bicu_stencil(1,14) =                 &
            get_gauss_opposite_neighbour_cell(                   &
               base(1), g_nbr_points_per_lat(global_num_lats),   &
                        g_nbr_points_per_lat(global_num_lats-1)) &
            + shift(1)

         psmile_gauss_shift_bicu_stencil(2,14) = global_num_lats-1

         orientation(4) = -1
   else if (base(2)+shift(2) == global_num_lats-1) then

         psmile_gauss_shift_bicu_stencil(1,14) =                 &
            get_gauss_opposite_neighbour_cell(                   &
               base(1), g_nbr_points_per_lat(global_num_lats-1), &
                        g_nbr_points_per_lat(global_num_lats))   &
            + shift(1)

         psmile_gauss_shift_bicu_stencil(2,14) = global_num_lats

         orientation(4) = -1
   else

         psmile_gauss_shift_bicu_stencil(1,14) =                               &
            get_gauss_neighbour_cell(base(1),                                  &
                                     g_nbr_points_per_lat(base(2)),            &
                                     g_nbr_points_per_lat(base(2)+2+shift(2))) &
            + shift(1)

         psmile_gauss_shift_bicu_stencil(2,14) = base(2) + 2 + shift(2)
   endif

   psmile_gauss_shift_bicu_stencil(3,14) = base(3) + shift(3)

   ! next we compute the longitude indices of the other columns

   do i = 2, 4

      ! same as i_ = mod(i,4)+1
      ! this will generate the indices 2,3,1
      i_ = iand (i, 3) + 1

      psmile_gauss_shift_bicu_stencil(1,i_::4) = psmile_gauss_shift_bicu_stencil(1,2::4) &
                                               + lon_shift(i-1) * orientation(:)
      psmile_gauss_shift_bicu_stencil(2:3,i_::4) = psmile_gauss_shift_bicu_stencil(2:3,2::4)
   enddo

   ! last we need to do some periodic index adjustment

   where (psmile_gauss_shift_bicu_stencil(1,:) < 1)
      psmile_gauss_shift_bicu_stencil(1,:) = psmile_gauss_shift_bicu_stencil(1,:) + &
                                           g_nbr_points_per_lat(psmile_gauss_shift_bicu_stencil(2,:))
   end where

   where (psmile_gauss_shift_bicu_stencil(1,:) > &
          g_nbr_points_per_lat(psmile_gauss_shift_bicu_stencil(2,:)))
      psmile_gauss_shift_bicu_stencil(1,:) = psmile_gauss_shift_bicu_stencil(1,:) - &
                                           g_nbr_points_per_lat(psmile_gauss_shift_bicu_stencil(2,:))
   end where

end function psmile_gauss_shift_bicu_stencil

!--------------------------------------------------------------------------------------------------

function psmile_gauss_get_neigh_stencil (grid_id, base)

   use psmile_common, only             : ndim_3d
   use psmile, only                    : grids
   use psmile_grid_reduced_gauss, only : get_gauss_opposite_neighbour_cell, &
                                         get_gauss_neighbour_cell

   implicit none

   integer, intent (in) :: grid_id, base(ndim_3d)

   integer :: psmile_gauss_get_neigh_stencil(ndim_3d, 9)

   integer            :: global_num_lats, i
   integer, pointer   :: g_nbr_points_per_lat(:)
   integer, parameter :: lon_shift(3) = (/-1,0,1/)

   ! the points returned have the following relation to each other ("base" has index 5):
   !
   !     7   8   9
   !
   !     4   5   6
   !
   !     1   2   3

   g_nbr_points_per_lat => grids(grid_id)%reduced_gauss_data%nbr_points_per_lat
   global_num_lats = size (g_nbr_points_per_lat)

   psmile_gauss_get_neigh_stencil(3, :) = base(3)
   psmile_gauss_get_neigh_stencil(:, 5) = base

   ! at first we compute the latitudes of all four rows

   psmile_gauss_get_neigh_stencil(2,1:3) = max (base(2)-1, 1)
      ! lat index of southern neighbours (if we are on the south pole
      ! then our neighbour is on the other side of the pole)
   psmile_gauss_get_neigh_stencil(2,4:6) = base(2)
      ! lat index of the base latitude
   psmile_gauss_get_neigh_stencil(2,7:9) = min (base(2)+1, global_num_lats)
      ! latitude index of northern neighbours (if we are on the
      ! north pole then our neighbour is on the other side
      ! of the pole)

   ! at second we compute the longitude indices of the base column (2,10,14)

   do i = 2, 8, 6

      ! if the point is supposed to be on the other side of the pole
      if (base(2) == psmile_gauss_get_neigh_stencil(2,i)) then

         psmile_gauss_get_neigh_stencil(1,i) =                           &
            get_gauss_opposite_neighbour_cell(base(1),                   &
                                          g_nbr_points_per_lat(base(2)), &
                                          g_nbr_points_per_lat(psmile_gauss_get_neigh_stencil(2,i)))
      else

         psmile_gauss_get_neigh_stencil(1,i) =                      &
            get_gauss_neighbour_cell(base(1),                       &
                                     g_nbr_points_per_lat(base(2)), &
                                     g_nbr_points_per_lat(psmile_gauss_get_neigh_stencil(2,i)))
      endif

   enddo

   ! next we compute the longitude indices of the other columns

   do i = 1, 3, 2

      psmile_gauss_get_neigh_stencil(1,i::3) = psmile_gauss_get_neigh_stencil(1,2::3) + lon_shift(i)
   enddo

   ! last we need to do some periodic index adjustment

   where (psmile_gauss_get_neigh_stencil(1,:) < 1)
      psmile_gauss_get_neigh_stencil(1,:) = psmile_gauss_get_neigh_stencil(1,:) + &
                                           g_nbr_points_per_lat(psmile_gauss_get_neigh_stencil(2,:))
   end where

   where (psmile_gauss_get_neigh_stencil(1,:) > &
          g_nbr_points_per_lat(psmile_gauss_get_neigh_stencil(2,:)))
      psmile_gauss_get_neigh_stencil(1,:) = psmile_gauss_get_neigh_stencil(1,:) - &
                                           g_nbr_points_per_lat(psmile_gauss_get_neigh_stencil(2,:))
   end where

end function psmile_gauss_get_neigh_stencil

!--------------------------------------------------------------------------------------------------

function psmile_gauss_get_bili_stencil (grid_id, base)

   use psmile_common, only             : ndim_3d
   use psmile_grid_reduced_gauss, only : psmile_gauss_shift_bili_stencil

   implicit none

   integer, intent (in) :: grid_id, base(ndim_3d)

   integer :: psmile_gauss_get_bili_stencil(ndim_3d, 4)


   psmile_gauss_get_bili_stencil = psmile_gauss_shift_bili_stencil (grid_id, base, (/0,0,0/))

end function psmile_gauss_get_bili_stencil

!--------------------------------------------------------------------------------------------------

function psmile_gauss_shift_bili_stencil (grid_id, base, shift)

   use psmile_common, only             : ndim_3d
   use psmile, only                    : grids, psmile_assert
   use psmile_grid_reduced_gauss, only : get_gauss_opposite_neighbour_cell, &
                                         get_gauss_neighbour_cell

   implicit none

   integer, intent (in) :: grid_id, base (ndim_3d), shift(ndim_3d)

   integer :: psmile_gauss_shift_bili_stencil (ndim_3d, 4)

   integer            :: global_num_lats
   integer, pointer   :: g_nbr_points_per_lat(:)
   integer            :: orientation

   ! the points returned have the following relation to each other ("base" has index 6):
   !
   !     4   3
   !
   !     1   2

   g_nbr_points_per_lat => grids(grid_id)%reduced_gauss_data%nbr_points_per_lat
   global_num_lats = size (g_nbr_points_per_lat)
   orientation = 1

#ifdef PRISM_ASSERTION
   if (any (abs (shift(:)) > 1)) then
      call psmile_assert (__FILE__, __LINE__, &
                          "shifting by more than one cell is not supported")
   endif
   if (base(2)+shift(2) < 1 .or. &
       base(2)+shift(2) > size (g_nbr_points_per_lat)) then
      call psmile_assert (__FILE__, __LINE__, &
                          "shifting of the base point over the pole is not supported")
   endif
#endif

   ! compute the first point

   ! if we need to recalculate the i index of point 1
   if (abs(shift(2)) == 1) then
   
      psmile_gauss_shift_bili_stencil(1,1) =                                         &
         get_gauss_neighbour_cell(base(1), g_nbr_points_per_lat(base(2)),          &
                                           g_nbr_points_per_lat(base(2)+shift(2))) &
         + shift(1)

      psmile_gauss_shift_bili_stencil(2,1) = base(2) + shift(2)

   else ! shift(2) == 0

      psmile_gauss_shift_bili_stencil(1:2,1) = base(1:2) + shift(1:2)
   endif

   psmile_gauss_shift_bili_stencil(3,1) = base(3) + shift(3)

   ! compute the fourth point

   ! if we need to recalculate the i index of point 4
   if (shift(2) >= 0) then

      ! some special care needs to be taken in case it is on the other side of the pole
      if (base(2)+shift(2) == global_num_lats) then

            psmile_gauss_shift_bili_stencil(1,4) =                &
               get_gauss_opposite_neighbour_cell(                 &
                  base(1), g_nbr_points_per_lat(global_num_lats), &
                           g_nbr_points_per_lat(global_num_lats)) &
               + shift(1)

            psmile_gauss_shift_bili_stencil(2,4) = global_num_lats

            orientation = -1
      else

            psmile_gauss_shift_bili_stencil(1,4) =                                &
               get_gauss_neighbour_cell(base(1),                                  &
                                        g_nbr_points_per_lat(base(2)),            &
                                        g_nbr_points_per_lat(base(2)+1+shift(2))) &
               + shift(1)

            psmile_gauss_shift_bili_stencil(2,4) = base(2) + 1 + shift(2)
      endif

   else ! shift(2) == -1
   
      psmile_gauss_shift_bili_stencil(1,4) = base(1) + shift(1)
      psmile_gauss_shift_bili_stencil(2,4) = base(2)
   endif

   psmile_gauss_shift_bili_stencil(3,4) = base(3) + shift(3)

   ! compute the second and third point
   psmile_gauss_shift_bili_stencil(:,2) = psmile_gauss_shift_bili_stencil(:,1) + (/1,0,0/) 
   psmile_gauss_shift_bili_stencil(:,3) = psmile_gauss_shift_bili_stencil(:,4) + &
                                          (/orientation,0,0/)
   
   ! last we need to do some periodic index adjustment
   where (psmile_gauss_shift_bili_stencil(1,:) < 1)
      psmile_gauss_shift_bili_stencil(1,:) = psmile_gauss_shift_bili_stencil(1,:) + &
                                          g_nbr_points_per_lat(psmile_gauss_shift_bili_stencil(2,:))
   end where

   where (psmile_gauss_shift_bili_stencil(1,:) > &
         g_nbr_points_per_lat(psmile_gauss_shift_bili_stencil(2,:)))
      psmile_gauss_shift_bili_stencil(1,:) = psmile_gauss_shift_bili_stencil(1,:) - &
                                          g_nbr_points_per_lat(psmile_gauss_shift_bili_stencil(2,:))
   end where

end function psmile_gauss_shift_bili_stencil

!--------------------------------------------------------------------------------------------------

function psmile_gauss_3d_to_global_1d (grid_id, points_3d, num_points)

   use psmile_common, only : ndim_3d
   use psmile, only        : grids

   implicit none

   integer, intent(in) :: grid_id, num_points
   integer, intent(in) :: points_3d(ndim_3d, num_points)

   integer :: psmile_gauss_3d_to_global_1d(num_points)

   integer :: i, num_lats

   num_lats = size (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat)

   do i = 1, num_points

      if ((points_3d(2,i) == 1) .or. (points_3d(2,i) > num_lats)) then

         psmile_gauss_3d_to_global_1d(i) = points_3d(1,i)

      else

         if (points_3d(1,i) >= 1 .and. &
             points_3d(1,i) <= grids(grid_id)%reduced_gauss_data%nbr_points_per_lat(points_3d(2,i))) then
             
            psmile_gauss_3d_to_global_1d(i) = points_3d(1,i) + &
               sum (grids(grid_id)%reduced_gauss_data%nbr_points_per_lat(1:points_3d(2,i)-1))
         else
            psmile_gauss_3d_to_global_1d(i) = points_3d(1,i)
         endif
      endif
   enddo

end function psmile_gauss_3d_to_global_1d

!--------------------------------------------------------------------------------------------------

function psmile_gauss_1d_global_to_local (grid_id, points_1d, num_points, &
                                          fill_value, inc_remote_neigh)

   use psmile, only        : grids

   implicit none

   integer, intent(in)           :: grid_id, num_points, fill_value
   integer, intent(in)           :: points_1d(num_points)
   logical, intent(in), optional :: inc_remote_neigh ! is this routine supposed to
                                                     ! include the remote points of
                                                     ! which this processes has a
                                                     ! local copy

   integer :: psmile_gauss_1d_global_to_local(num_points)

   integer :: i, j, num_local_points, curr_block

   psmile_gauss_1d_global_to_local(:) = fill_value

   if (associated (grids(grid_id)%reduced_gauss_data%local_block_info)) then

      curr_block = 1

      ! for all points
   outer: do i = 1, num_points

         do j = curr_block, size (grids(grid_id)%partition, 1)

            if (points_1d(i) >= &
                grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_start_idx &
                .and. &
                points_1d(i) <= &
                grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_end_idx) then

               psmile_gauss_1d_global_to_local(i) = &
                  grids(grid_id)%reduced_gauss_data%local_block_info(j)%local_1d_start_idx + &
                  points_1d(i) - &
                  grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_start_idx

               curr_block = j
               cycle outer
            endif
         enddo

         do j = 1, curr_block-1

            if (points_1d(i) >= &
                grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_start_idx &
                .and. &
                points_1d(i) <= &
                grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_end_idx) then

               psmile_gauss_1d_global_to_local(i) = &
                  grids(grid_id)%reduced_gauss_data%local_block_info(j)%local_1d_start_idx + &
                  points_1d(i) - &
                  grids(grid_id)%reduced_gauss_data%local_block_info(j)%global_1d_start_idx

               curr_block = j
               cycle outer
            endif
         enddo
      enddo outer ! i

   else ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

      ! for all blocks in the local partition
      do i = 1, size (grids(grid_id)%partition, 1)
      
         ! look for points that are in the current block and transform them to local indices
         where (points_1d(:) > grids(grid_id)%partition(i,1) .and. &
               points_1d(:) <= grids(grid_id)%partition(i,1) + grids(grid_id)%extent(i,1))
      
            psmile_gauss_1d_global_to_local(:) = points_1d(:) - &
                                                grids(grid_id)%partition(i,1) + &
                                                sum (grids(grid_id)%extent(1:i-1,1))
      
         end where
      
      enddo

   endif ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

   ! if the user specified that the locally stored remote points
   ! are not to be taken into account -> return
   if (present(inc_remote_neigh)) then
      if (.not. inc_remote_neigh) return
   endif

   ! if we have some data on remote points stored locally
   if (associated (grids(grid_id)%remote_index)) then

      num_local_points = sum (grids(grid_id)%extent(:,1))

      ! for all points which were not found in the local partitions
outer2:do i = 1, num_points

         if (psmile_gauss_1d_global_to_local(i) == fill_value) then

            do j = 1, size (grids(grid_id)%remote_index)

               if (points_1d(i) == grids(grid_id)%remote_index(j)) then

                  psmile_gauss_1d_global_to_local(i) = j + num_local_points
                  cycle outer2

               endif
            enddo ! j
         endif
      enddo outer2 ! i
   endif

end function psmile_gauss_1d_global_to_local

!--------------------------------------------------------------------------------------------------

function psmile_gauss_3d_to_local_1d (grid_id, points_3d, num_points, &
                                          fill_value, inc_remote_neigh)

   use psmile, only                    : grids
   use psmile_common, only             : ndim_3d
   use psmile_grid_reduced_gauss, only : psmile_gauss_3d_to_global_1d, &
                                         psmile_gauss_1d_global_to_local, &
                                         block_info

   implicit none

   integer, intent(in)           :: grid_id, num_points, fill_value
   integer, intent(in)           :: points_3d(ndim_3d, num_points)
   logical, intent(in), optional :: inc_remote_neigh ! is this routine supposed to
                                                     ! include the remote points of
                                                     ! which this processes has a
                                                     ! local copy

   integer :: psmile_gauss_3d_to_local_1d(num_points)

   integer                    :: i, j, num_blocks, num_local_points
   integer                    :: point_global_1d(1)
   type (block_info), pointer :: local_block_info (:)

   ! if there is additional information on the local blocks available we can
   ! optimise the transformation
   if (associated (grids(grid_id)%reduced_gauss_data%local_block_info)) then

       psmile_gauss_3d_to_local_1d = fill_value
       local_block_info => grids(grid_id)%reduced_gauss_data%local_block_info
       num_blocks = size (local_block_info)

!       do n = 1, num_points
!         do i = 1 , num_blocks
 
!            if (points_3d(1, n) >= local_block_info(i)%global_3d_start_idx(1) .and. &
!                points_3d(1, n) <= local_block_info(i)%global_3d_end_idx(1) .and.   &
!                points_3d(2, n) == local_block_info(i)%global_3d_start_idx(2)) then
 
!               psmile_gauss_3d_to_local_1d(n) =                &
!                  local_block_info(i)%local_1d_start_idx +     &
!                  points_3d(1, n) -                            &
!                  local_block_info(i)%global_3d_start_idx(1)
!               exit
!            endif    
!         enddo ! i
!       enddo ! n
      do i = 1 , num_blocks

         where (points_3d(1, :) >= local_block_info(i)%global_3d_start_idx(1) .and. &
                points_3d(1, :) <= local_block_info(i)%global_3d_end_idx(1) .and.   &
                points_3d(2, :) == local_block_info(i)%global_3d_start_idx(2))

            psmile_gauss_3d_to_local_1d(:) =                &
               local_block_info(i)%local_1d_start_idx +     &
               points_3d(1, :) -                            &
               local_block_info(i)%global_3d_start_idx(1)
         endwhere    
      enddo ! i

      if (present(inc_remote_neigh)) then

         if (inc_remote_neigh) then

            num_local_points = sum (grids(grid_id)%extent(:,1))

outer:      do i = 1, num_points

               if  (psmile_gauss_3d_to_local_1d(i) == fill_value) then

                  point_global_1d = psmile_gauss_3d_to_global_1d(grid_id, (/points_3d(:, i)/), 1)

                  do j = 1, size (grids(grid_id)%remote_index)

                     if (point_global_1d(1) == grids(grid_id)%remote_index(j)) then

                        psmile_gauss_3d_to_local_1d(i) = j + num_local_points
                        cycle outer

                     endif
                  enddo ! j
               endif
            enddo outer
         endif ! inc_remote_neigh
      endif ! present(inc_remote_neigh)

   else ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

      if (present(inc_remote_neigh)) then
         psmile_gauss_3d_to_local_1d = psmile_gauss_1d_global_to_local (grid_id,     &
                                          psmile_gauss_3d_to_global_1d (grid_id,     &
                                                                        points_3d,   &
                                                                        num_points), &
                                          num_points, fill_value, inc_remote_neigh)
      else
         psmile_gauss_3d_to_local_1d = psmile_gauss_1d_global_to_local (grid_id,     &
                                          psmile_gauss_3d_to_global_1d (grid_id,     &
                                                                        points_3d,   &
                                                                        num_points), &
                                          num_points, fill_value)
      endif
   endif ! associated (grids(grid_id)%reduced_gauss_data%local_block_info)

end function psmile_gauss_3d_to_local_1d