psmile_mg_method_2d_dble.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! Copyright 2010, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_MG_method_2d_dble
!
! !INTERFACE:

subroutine psmile_mg_method_2d_dble (                    &
                comp_info, nlev, found, loc, range,      &
                coords1, coords2, search_shape, control, &
                x_coords, y_coords, coords_shape,        &
                grid_valid_shape,                        &
                cyclic, period,                          &
                chmin1, chmin2, chmax1, chmax2,          &
                tol, ierror)
!
! !USES:
!
  use prism_constants
  use psmile, dummy_interface => PSMILe_mg_method_2d_dble
  use psmile_geometry, only : point_is_in_quadrangle
  use psmile_grid, only : pole_threshold, psmile_transform_cell_cyclic
#ifdef DEBUG_TRACE
  use psmile_debug_trace
#endif

  implicit none

!
! !INPUT PARAMETERS:
!

  Type (Enddef_comp), Intent (In) :: comp_info

! Info on the component in which the donor should be searched.

  Integer, Intent (In)            :: nlev

! Total number of levels

  Integer, Intent (In)            :: range (2, ndim_3d)

! Dimension of loc and found

  Integer, Intent (In)            :: search_shape (2, ndim_3d)

! Dimension of coordinate arrays "coord1" and "coord2"
! which contain the coordinates to be searched.

  Double Precision, Intent (In)   :: coords1(search_shape(1,1):search_shape(2,1), 
                                             search_shape(1,2):search_shape(2,2), 
                                             search_shape(1,3):search_shape(2,3))

  Double Precision, Intent (In)   :: coords2(search_shape(1,1):search_shape(2,1), 
                                             search_shape(1,2):search_shape(2,2), 
                                             search_shape(1,3):search_shape(2,3))
! Coordinates to be searched

  Integer, Intent (In)            :: control (2, ndim_3d)

! Index range in "coords" to be searched

  Integer, Intent (In)            :: coords_shape (2, ndim_2d)

! Dimension of coordinates (method) x_coords, ...

  Double Precision, Intent (In)   :: x_coords(coords_shape(1,1):coords_shape(2,1), 
                                              coords_shape(1,2):coords_shape(2,2))

  Double Precision, Intent (In)   :: y_coords(coords_shape(1,1):coords_shape(2,1), 
                                              coords_shape(1,2):coords_shape(2,2))

! Coordinates of the method

  Integer, Intent (In)            :: grid_valid_shape (2, ndim_2d)

! Specifies the valid block shape for the "ndim_2d"-dimensional block
! Dimension of chmin, chmax and midp

  Logical,          Intent (In)   :: cyclic (ndim_2d)

! Is this a (locally) cyclic coordinate in corresponding direction

  Double Precision, Intent (In)   :: period (ndim_2d)
!
! Period of a cyclic coordinate.
!
  Double Precision, Intent (In)   :: chmin1 (grid_valid_shape(1,1)-1:grid_valid_shape(2,1), 
                                             grid_valid_shape(1,2)-1:grid_valid_shape(2,2))

  Double Precision, Intent (In)   :: chmin2 (grid_valid_shape(1,1)-1:grid_valid_shape(2,1), 
                                             grid_valid_shape(1,2)-1:grid_valid_shape(2,2))

! Minimum of grid coordinates per grid cell
! chmin (i, j) contains the minimal value of bounding box of
!  method cell ((i,j), (i+1,j), (i,j+1), (i+1,j+1))
!
! The first box
!     chmin (grid_valid_shape(1,1)-1,j) contains the minimal value of
!                                   the bounding box of virtual cell
!       (array(grid_valid_shape(1,1), j), array(grid_valid_shape(1,1), j+1),
!        bounding box of corner control volume of (grid_valid_shape(1,1),j)
!
! The last box
!     chmin (grid_valid_shape(2,1),  *) contains the minimal value of
!                                   the bounding box of virtual cell
!       (array(grid_valid_shape(2,1), j), array(grid_valid_shape(2,1), j+1),
!        bounding box of corner control volume of (grid_valid_shape(2,1),j)
!
! Corresponding values are contained for other virtual cells at block
! border.
!
  Double Precision, Intent (In)   :: chmax1 (grid_valid_shape(1,1)-1:grid_valid_shape(2,1), 
                                             grid_valid_shape(1,2)-1:grid_valid_shape(2,2))
  Double Precision, Intent (In)   :: chmax2 (grid_valid_shape(1,1)-1:grid_valid_shape(2,1), 
                                             grid_valid_shape(1,2)-1:grid_valid_shape(2,2))


! Maximum of grid coordinates per grid cell

  Double Precision, Intent (In)   :: tol
!
! Absolute tolerance for search of "identical" points TOL >= 0.0
!
! !INPUT/OUTPUT PARAMETERS:
!
  Integer, Intent (InOut)         :: found (range(1,1):range(2,1), 
                                            range(1,2):range(2,2), 
                                            range(1,3):range(2,3))

! Finest level number on which a grid cell was found for point i,j,k.
! Input :
!     Level number < -1: Point was not found and
!                        and last level number was (-found(i,j,k))
!     Level number = -(nlev+1): Never found
! Output :
!     found(i,j,k) = +1: Point (i,j,k) is located on a point
!                        of the method grid.
!     found(i,j,k) = -1: Point (i,j,k) is located in a cell
!                        of the method grid.
!     Otherwise, not contained in the method grid.
!
  Integer, Intent (InOut)         :: loc (ndim_2d,               
                                          range(1,1):range(2,1), 
                                          range(1,2):range(2,2), 
                                          range(1,3):range(2,3))

! Indices of the grid cell in which the point was found.
! The indices are relative to "shape".
!
! !OUTPUT PARAMETERS:
!
  Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_mg_method_2d_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  val_direct   = Code for locations which should to be directly transferred
!                 to the destination process.
!  val_coupler  = Code for locations which should to be transferred
!                 to the coupler process.
!  val_notfound = Code for locations which was not found
!
  Integer, Parameter              :: nrefd = 3 * 3
  Integer, Parameter              :: n_corners = 4

  Integer, Parameter              :: val_direct   =  1
  Integer, Parameter              :: val_coupler  = -1
  Integer                         :: val_notfound

  Integer, Parameter              :: temp_len = 256
!
!     ... for levels
!
  Integer, Parameter              :: lev = 1
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
  Integer                         :: i, j, k, l
  Integer                         :: ibegl, jpart, n
  Integer                         :: curr_loc_list (temp_len, ndim_2d)
!
!     ... For locations controlled
!
  Integer                         :: ijkctl (ndim_2d, nrefd)
  Integer                         :: ijkdst (temp_len, ndim_2d, nrefd)
  Integer                         :: irange (temp_len, 2, ndim_2d)

  Integer                         :: ii, jj, nref (temp_len)
!
!     ... For cyclic cells
!
  Logical                         :: cyclic_req

  Double Precision                :: period2 (ndim_2d)
!
!     ... Local reference point and edges
!
  Type (point_dble)               :: xy_source(temp_len, 4)
  Type (point_dble)               :: xy_source_point(temp_len)
  Type (point_dble), target       :: xy_target(temp_len)
  Type (point_dble), target       :: xy_target_rotated(temp_len)
  Type p_xy_target
     Type (point_dble), pointer   :: p
  end Type p_xy_target
  Type(p_xy_target)               :: xy_target_pointer(temp_len)
!
!     ...
!
  Integer                         :: num_curr_locs, nc_full, np, nc_part
  Integer                         :: list (temp_len*2), full_list (temp_len)
  Integer                         :: part_list (temp_len)
  Logical                         :: fnd (temp_len)

#ifdef DEBUG
!     ... for locations searched
!
  Integer                         :: nfnd (0:3)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_method_2d_dble" searches the corresponding points
! of the method (grid) for the subgrid coords by the sending process.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
! 10.09.10    M. Hanke     major revision to address problem at poles, date line, and 0-Meridian
!             R. Redler
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_method_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $
!  $Author: hanke $
!
  Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_method_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices
!
  data ((ijkctl (i, n), i=1,ndim_2d), n = 1,nrefd) &
       / -1,-1,   0,-1,   1,-1, &
       -1, 0,   0, 0,   1, 0, &
       -1, 1,   0, 1,   1, 1/
!----------------------------------------------------------------------
!
!  Initialization
!
  ierror = 0
  val_notfound = -(nlev+1)

#ifdef VERBOSE
  print 9990, trim(ch_id), lev, control

  call psmile_flushstd
#endif /* VERBOSE */

#ifdef DEBUG_TRACE
  if ( range(1,1) <= ictl_ind(1) .and. range(2,1) >= ictl_ind(1) .and. &
       range(1,2) <= ictl_ind(2) .and. range(2,2) >= ictl_ind(2)) then
     print 8990, trim(ch_id), ictl_ind (1:ndim_2d), &
          found ( ictl_ind(1),ictl_ind(2),control(1,3)), &
          loc (:, ictl_ind(1),ictl_ind(2),control(1,3))
  endif
#endif /* DEBUG_TRACE */
  !
#ifdef PRISM_ASSERTION
  if (tol < 0.0) then
     call psmile_assert (__FILE__, __LINE__, &
          "tol must be >= 0.0")
  endif

  if ( grid_valid_shape(1,1)-1 < coords_shape (1,1) .or. &
       grid_valid_shape(2,1)+1 > coords_shape (2,1) .or. &
       grid_valid_shape(1,2)-1 < coords_shape (1,2) .or. &
       grid_valid_shape(2,2)+1 > coords_shape (2,2)) then

     print *, 'grid_valid_shape:', grid_valid_shape
     print *, 'coords_shape    :', coords_shape
     call psmile_assert (__FILE__, __LINE__, &
          "insufficient coords_shape")
  endif
#endif
  !
  ! create temporary memory (control (2,1)-control(1,1)+1)*ndtry
  ! and do n =1,nref before i-loop ?!?
  !
  cyclic_req = cyclic (1) .or. cyclic (2)
  period2 (:) = period (:) * 0.5
  !
#ifdef DEBUG
  nfnd (:) = 0
#endif /* DEBUG */
  !
  !===>
  !

  do k = control(1,3), control(2,3)
     do j = control(1,2), control (2,2)

!
!-----------------------------------------------------------------------
!             Look in the fine parts of the coarse grid cell
!             loc(:, i,j,k)
!
!            TODO: List over 1-dimensional Vector (i,j,k)
!-----------------------------------------------------------------------
!
        ibegl = control (1, 1)
!
!===> ... Look for the next cells to be controlled
!         Store indices of cells in list "list"
!
loop_i: do while (ibegl <= control(2,1))

           num_curr_locs = 0
           i = ibegl

           do while (num_curr_locs < temp_len .and. i <= control(2,1))
              ibegl = i
!cdir vector
              do i = ibegl, min (ibegl+(temp_len*2-1)-num_curr_locs, control(2,1))
                 if (found(i, j, k) == lev) then
                    num_curr_locs = num_curr_locs + 1
                    list (num_curr_locs) = i
                 endif
              end do
           end do

           if (num_curr_locs == 0) exit loop_i ! exit while loop (ibegl < ...)

           if (num_curr_locs < temp_len) then
              ibegl = control(2,1) + 1
           else  if (num_curr_locs > temp_len) then
              num_curr_locs = temp_len
              ibegl = list (temp_len+1)
           else
              ibegl = list (temp_len) + 1
           endif

           if (cyclic_req) then
              nc_full = num_curr_locs
              full_list (1:num_curr_locs) = list (1:num_curr_locs)
           else
              nc_full = 0
           endif

!
!===> ... Look in the real cell (curr_loc_list(i, 1),curr_loc_list(i, 2))
!
! get the current locations
           curr_loc_list (1:num_curr_locs, 1) = loc (1, list(1:num_curr_locs), j, k)
           curr_loc_list (1:num_curr_locs, 2) = loc (2, list(1:num_curr_locs), j, k)

#ifdef PRISM_ASSERTION
!cdir vector
           do i = 1, num_curr_locs
              if ( curr_loc_list(i,1) < coords_shape(1,1) .or. &
                   curr_loc_list(i,1) > coords_shape(2,1) .or. &
                   curr_loc_list(i,2) < coords_shape(1,2) .or. &
                   curr_loc_list(i,2) > coords_shape(2,2) ) exit
           end do

           if (i <= num_curr_locs) then
              print *, 'i, curr_loc_list', i, curr_loc_list(i, :)
              print *, 'coords_shape', coords_shape
              call psmile_assert (__FILE__, __LINE__, &
                   'wrong index')
           endif
#endif
!cdir vector
           ! for each location in the current list
           do i = 1, num_curr_locs
              ! get the range which contains the neighbours of the current location.
              ! the range is trunctated at the grid boundaries
              irange (i, 1, :) = max (curr_loc_list(i,:)-1, grid_valid_shape (1, :))
              irange (i, 2, :) = min (curr_loc_list(i,:)+1, grid_valid_shape (2, :))
           end do

           ! compute the number of available neighbours + the actual location
           nref (1:num_curr_locs) = &
               (irange (1:num_curr_locs, 2,1) - irange (1:num_curr_locs, 1,1) + 1) * &
               (irange (1:num_curr_locs, 2,2) - irange (1:num_curr_locs, 1,2) + 1)

#ifdef PRISM_ASSERTION
           if (minval (nref(1:num_curr_locs)) <= 0) then
              call psmile_assert (__FILE__, __LINE__, &
                   'nref <= 0')
           endif
#endif
!cdir vector
           do i = 1, num_curr_locs
              ! if the all (9) neighbours were available (no trunctation at grid boundary)
              if (nref(i) == nrefd) then
                 do n = 1, nrefd
                    ! generate list of all neighbours + current location
                    ijkdst (i, :, n) = curr_loc_list (i, :) + ijkctl (:, n)
                 end do
              else
                 ! generate list of all neighbours + current location
                 n = 0
                 do jj = irange (i, 1, 2), irange (i, 2, 2)
                    do ii = irange (i, 1, 1), irange (i, 2, 1)
                       n = n + 1
                       ijkdst (i, 1, n) = ii
                       ijkdst (i, 2, n) = jj
                    end do ! ii
                 end do ! jj
              endif
           end do

           ! get the target coords for all locations currently processed
!cdir vector
           do i = 1, num_curr_locs
              xy_target(i)%x = coords1(list(i),j,k)
              xy_target(i)%y = coords2(list(i),j,k)
           enddo

           nc_part = 0

!
!===> ... Check for matching points
!
           do n = 1, nrefd
!cdir vector
              do i = 1, num_curr_locs
                 if (n <= nref(i)) then
                    xy_source_point%x = x_coords (ijkdst(i,1,n), ijkdst (i,2,n))
                    xy_source_point%y = y_coords (ijkdst(i,1,n), ijkdst (i,2,n))
                    ! if the source cell is close to the pole
                    ! transform source and target point to the equator
                    if ( y_coords (ijkdst(i,1,n), ijkdst (i,2,n)) >  pole_threshold .or. &
                         y_coords (ijkdst(i,1,n), ijkdst (i,2,n)) < -pole_threshold ) then

                       call psmile_transrot_dble (                     &
                             x_coords (ijkdst(i,1,n), ijkdst (i,2,n)), &
                             y_coords (ijkdst(i,1,n), ijkdst (i,2,n)), &
                             xy_source_point(i)%x, xy_source_point(i)%y )

                       xy_target_rotated(i) = xy_target(i)

                       call psmile_transrot_dble (                        &
                             coords1(list(i),j,k), coords2(list(i),j,k),   &
                             xy_target_rotated(i)%x, xy_target_rotated(i)%y )

                       xy_target_pointer(i)%p => xy_target_rotated(i)
                    else
                       xy_source_point(i)%x = x_coords (ijkdst(i,1,n), ijkdst (i,2,n))
                       xy_source_point(i)%y = y_coords (ijkdst(i,1,n), ijkdst (i,2,n))
                       xy_target_pointer(i)%p => xy_target(i)
                    endif

                    ! if source and target point match
                    if (abs(xy_target_pointer(i)%p%x - xy_source_point(i)%x) < tol .and. &
                        abs(xy_target_pointer(i)%p%y - xy_source_point(i)%y) < tol) then
!not needed            found (list(i),j,k)  = val_direct ! == lev == 1
                       loc (:, list(i),j,k) = ijkdst (i, :, n)
                       nc_part = nc_part + 1
                       part_list (nc_part) = i
                    endif
                 end if
              end do
           end do

           ! if matching points were found
           if (nc_part > 0) then

#ifdef DEBUG
              nfnd (1) = nfnd (1) + nc_part
#endif /* DEBUG */

              ! if all target cells in the current list have a matching source cell
              ! get the next list of target points
              if (num_curr_locs == nc_part) cycle loop_i
!
!===> ...... Remove indices of "part_list" from "list"
!            ??? Vectorisierung mittels hilfsvektors ?
!
              do np = 0, nc_part-1
                 i = part_list (nc_part - np)

                 if (i /= num_curr_locs-np) then
                    list (i) = list (num_curr_locs-np)
                    nref (i) = nref (num_curr_locs-np)
                    xy_target(i) = xy_target(num_curr_locs-np)

                    curr_loc_list (i, :) = curr_loc_list (num_curr_locs-np, :)

                    ijkdst (i, :, :) = ijkdst (num_curr_locs-np, :, :)
                    irange (i, :, :) = irange (num_curr_locs-np, :, :)
                 endif
              end do

              num_curr_locs = num_curr_locs - nc_part
           endif
!
!-------------------------------------------------------------------------------
!     ... Points are not located on the grid point of the method.
!         Look for the correct method cell.
!-------------------------------------------------------------------------------
!
!cdir vector
           do i = 1, num_curr_locs
              found (list(i), j, k) = val_coupler
           end do
!
!===> ....... Get cell range "irange" for search of points
!             which are located at the lower boundary
!             (irange (i,1,:) must be extended) (see virtual cells)
!
           nc_part = 0
!cdir vector
           do i = 1, num_curr_locs
              if ( curr_loc_list(i,1) == grid_valid_shape(1,1) .or. &
                   curr_loc_list(i,2) == grid_valid_shape(1,2)) then
                 nc_part = nc_part + 1
                 part_list (nc_part) = i
              endif
           end do
           !
!cdir vector
           do np = 1, nc_part
              i = part_list (np)
              !
              ! This will place some locations out of the grid_valid_shape. This must therefore be
              !  considered during the ASSERTION check in psmile_neigh_trili_irreg2
              !                                           psmile_neigh_near_irreg2
              !                                           psmile_neigh_tricu_irreg2, etc.
              !
              irange (i,1, :) = max (curr_loc_list(i,:)-1, grid_valid_shape(1,:)-1)
              irange (i,2, :) = min (curr_loc_list(i,:)+1, grid_valid_shape(2,:))

              nref(i) = (irange (i,2,1) - irange (i,1,1) + 1) * &
                        (irange (i,2,2) - irange (i,1,2) + 1)

              if (nref(i) == nrefd) then
                 do n = 1, nrefd
                    ijkdst (i, :, n) = curr_loc_list (i, :) + ijkctl (:, n)
                 end do
              else
                 n = 0
                 do jj = irange (i,1, 2), irange (i,2, 2)
                    do ii = irange (i,1, 1), irange (i,2, 1)
                       n = n + 1
                       ijkdst (i, 1, n) = ii
                       ijkdst (i, 2, n) = jj
                    end do ! ii
                 end do ! jj
              endif
           end do

#ifdef PRISM_ASSERTION
!cdir vector
           do np = 1, nc_part
              if (nref(part_list (np)) <= 0) exit
           end do

           if (np <= nc_part) then
              call psmile_assert (__FILE__, __LINE__, &
                   'nref <= 0')
           endif
#endif
!
!===> ... Find correct source method cell
!
           fnd (1:num_curr_locs) = .false.

           do n = 1, nrefd
!cdir vector
              do i = 1, num_curr_locs

                 ! This is the dead of vectorisation!
                 ! split loop over i into several smaller loops.
                 ! MoHa: once this runs we will have a look at
                 ! some performance measurments and determine
                 ! whether this is required or not

                 if ( n <= nref (i) .and. .not. fnd(i) ) then

                    ! move target longitudes into the common local longitude range
                    do while ( xy_target(i)%x < chmax1(ijkdst(i, 1,n),ijkdst(i, 2,n)) - period(1) )
                       xy_target(i)%x =  xy_target(i)%x + period(1)
                    enddo

                    do while ( xy_target(i)%x > chmin1(ijkdst(i, 1,n),ijkdst(i, 2,n)) + period(1) )
                       xy_target(i)%x =  xy_target(i)%x - period(1)
                    enddo

                    ! if the target coord is within the bounding box of the current cell
                    ! -> get the respective cell
                    if ( xy_target(i)%x >= chmin1(ijkdst(i,1,n), ijkdst(i,2,n)) .and. &
                         xy_target(i)%y >= chmin2(ijkdst(i,1,n), ijkdst(i,2,n)) .and. &
                         xy_target(i)%x <= chmax1(ijkdst(i,1,n), ijkdst(i,2,n)) .and. &
                         xy_target(i)%y <= chmax2(ijkdst(i,1,n), ijkdst(i,2,n))) then

#ifdef DEBUG_TRACE
                       if (j == ictl_ind(2)) then
                          if (list(i) == ictl_ind(1)) then
                             print *, 'is within bounding box:'
                             print *, 'ch1', chmin1(ijkdst(i,1,n), ijkdst(i,2,n)), &
                                               chmax1(ijkdst(i,1,n), ijkdst(i,2,n))
                             print *, 'ch2', chmin2(ijkdst(i,1,n), ijkdst(i,2,n)), &
                                               chmax2(ijkdst(i,1,n), ijkdst(i,2,n))
                          endif
                       end if
#endif /* DEBUG_TRACE */

                       ! if the bounding box is close to the pole
                       if ( chmin2(ijkdst(i,1,n),ijkdst(i,2,n)) >  pole_threshold .or. &
                            chmax2(ijkdst(i,1,n),ijkdst(i,2,n)) < -pole_threshold ) then

                          ! translate source and target grid coordintates to the equator
                          ! to do: make pole_threshold dependent on hor. resolution
                          call psmile_transrot_dble (                       &
                               x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  ), &
                               y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  ), &
                               xy_source(i, 1)%x, xy_source(i, 1)%y )

                          call psmile_transrot_dble (                       &
                               x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  ), &
                               y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  ), &
                               xy_source(i, 2)%x, xy_source(i, 2)%y )

                          call psmile_transrot_dble (                       &
                               x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1), &
                               y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1), &
                               xy_source(i, 3)%x, xy_source(i, 3)%y )

                          call psmile_transrot_dble (                       &
                               x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1), &
                               y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1), &
                               xy_source(i, 4)%x, xy_source(i, 4)%y )

                          call psmile_transrot_dble (                        &
                               coords1(list(i),j,k), coords2(list(i),j,k),   &
                               xy_target_rotated(i)%x, xy_target_rotated(i)%y )

                          ! the transformation might cause some bad coordinates
                          call psmile_transform_cell_cyclic ( xy_source(i, :)%x, period(1), ierror )

                          ! move target longitudes into the common local longitude range
                          do while ( xy_target_rotated(i)%x < minval(xy_source(i, :)%x) - period2(1) )
                             xy_target_rotated(i)%x =  xy_target_rotated(i)%x + period(1)
                          enddo

                          do while ( xy_target_rotated(i)%x > maxval(xy_source(i, :)%x) + period2(1) )
                             xy_target_rotated(i)%x =  xy_target_rotated(i)%x - period(1)
                          enddo

                          xy_target_pointer(i)%p => xy_target_rotated(i)

#ifdef DEBUG_TRACE
                          if (j == ictl_ind(2)) then
                             if (list(i) == ictl_ind(1)) then

                                print *, 'point and cell are rotated to equator'
                                print "('  point coordinates before: ', 2f8.2)", xy_target(i)
                                print "('  cell  coordinates before: ', 8f8.2)", &
                                   x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  ),  &
                                   y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  ),  &
                                   x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  ),  &
                                   y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  ),  &
                                   x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1),  &
                                   y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1),  &
                                   x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1),  &
                                   y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1)
                             endif
                          end if
#endif /* DEBUG_TRACE */
                       else

                          xy_source(i, 1)%x = x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  )
                          xy_source(i, 1)%y = y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)  )

                          xy_source(i, 2)%x = x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  )
                          xy_source(i, 2)%y = y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)  )

                          xy_source(i, 3)%x = x_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1)
                          xy_source(i, 3)%y = y_coords (ijkdst(i,1,n)+1, ijkdst(i,2,n)+1)

                          xy_source(i, 4)%x = x_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1)
                          xy_source(i, 4)%y = y_coords (ijkdst(i,1,n)  , ijkdst(i,2,n)+1)

                          xy_target_pointer(i)%p => xy_target(i)

                          ! if any corner of the current cell is outside its own bounding box, then
                          ! the cell requires transforming
                          if (any(xy_source(i, :)%x < chmin1(ijkdst(i,1,n),ijkdst(i,2,n))) .or. &
                              any(xy_source(i, :)%x > chmax1(ijkdst(i,1,n),ijkdst(i,2,n)))) then
                             call psmile_transform_cell_cyclic ( xy_source(i, :)%x, period(1), ierror )
                          endif

                       endif

                       ! check whether the point is within the cell
                       if (point_is_in_quadrangle(xy_target_pointer(i)%p, xy_source(i,:))) then
                          loc(1,list(i),j,k) = ijkdst(i,1,n)
                          loc(2,list(i),j,k) = ijkdst(i,2,n)
                          fnd (i) = .true.
                       endif
#ifdef DEBUG_TRACE
                       if (j == ictl_ind(2)) then
                          if (list(i) == ictl_ind(1)) then

                             print "('  point coordinates: ', 2f8.2)", xy_target_pointer(i)%p
                             print "('  cell  coordinates: ', 8f8.2)", xy_source(i,:)

                             if (fnd (i)) then
                                print *, 'point is within this cell'
                             else
                                print *, 'point is not within this cell'
                             endif
                          endif
                       end if
                    else
                       if (j == ictl_ind(2)) then
                          if (list(i) == ictl_ind(1)) then
                             print *, 'is not within bounding box:'
                             print *, 'ch1', chmin1(ijkdst(i,1,n), ijkdst(i,2,n)), &
                                               chmax1(ijkdst(i,1,n), ijkdst(i,2,n))
                             print *, 'ch2', chmin2(ijkdst(i,1,n), ijkdst(i,2,n)), &
                                               chmax2(ijkdst(i,1,n), ijkdst(i,2,n))
                          endif
                       end if
#endif /* DEBUG_TRACE */
                    endif
                 endif
              end do ! i
           end do ! n

           ! check if there are point which did not find a matching cell
           ! mark them as not found
           do n = 1, nrefd
!cdir vector
              do i = 1, num_curr_locs
                 if ( n <= nref (i) .and. .not. fnd(i) ) then
                    found(list(i), j, k) = val_notfound
                 endif
              enddo
           enddo

#ifdef DEBUG_TRACE
           if (j == ictl_ind(2)) then
              do i = 1, num_curr_locs
                 if (list(i) == ictl_ind(1)) exit
              end do
              !
              if (i <= num_curr_locs) then
                 print *, 'mg_method: ictl_ind', ictl_ind(1:ndim_2d), &
                      ': ', fnd(i), coords1(list(i),j,k), coords2(list(i),j,k)
                 print *, ' .... found, i, j: ', fnd(i)
                 if ( fnd(i) ) then
                   print *, ' .... found, i, j: ', loc(1:2,list(i),j,k)
                 else
                   print *, ' did not find new i, j!'
                 endif
              endif
           end if
#endif /* DEBUG_TRACE */
!
!-----------------------------------------------------------------------------
!     ... All points supposed to be sent to the coupler are examined
!     ... Is there an periodic index
!-----------------------------------------------------------------------------
!

           do l = 1, ndim_2d
              if (cyclic(l)) then
                 !
                 ! This is an cyclic/periodic index;
                 ! Move points, which are located in virtual cell 0
                 ! into cell "grid_valid_shape(2,l)"
                 !
                 do jpart = 1, num_curr_locs, 5000
!cdir vector loopcnt=5000
                    do i = jpart, min (num_curr_locs, jpart+5000-1)
                       if (loc (l, list(i),j,k) < grid_valid_shape (1,l) &
                           .and. fnd(i)) then
                           loc (l, list(i),j,k) = grid_valid_shape (2,l)
                       endif
                    enddo

                 enddo ! jpart
              endif
           end do ! l
!
!===> ... Get next list of points for first logical index
!
        end do loop_i ! while
     end do ! j
  end do ! k

!
!===> All done
!
#ifdef DEBUG_TRACE
  if ( range(1,1) <= ictl_ind(1) .and. range(2,1) >= ictl_ind(1) .and. &
       range(1,2) <= ictl_ind(2) .and. range(2,2) >= ictl_ind(2)) then

     print *, 'coords', coords1(ictl_ind(1),ictl_ind(2),control(1,3)), &
          coords2(ictl_ind(1),ictl_ind(2),control(1,3))

     print 8990,  trim(ch_id), ictl_ind(1:ndim_2d), &
          found ( ictl_ind(1),ictl_ind(2),control(1,3)), &
          loc (:, ictl_ind(1),ictl_ind(2),control(1,3))
  endif

8990 format (1x, a, ': psmile_mg_method_2d_dble: ictl_ind', 2i5, &
          '; found', i3, ' loc ', 2i6)
#endif /* DEBUG_TRACE */

#ifdef DEBUG
  print 9970, trim(ch_id), lev, nfnd
9970 format (1x, a, ': psmile_mg_method_2d_dble: lev =', i3, &
       ': not :', i5, ', dir: ', i6, ', fnd :', i6, i5)

#ifdef DEBUG_PRINT
  call psmile_print_found_dble (comp_info, lev,           &
       found, loc, range,                          &
       coords1, coords2, coords1,                  &
       search_shape, search_shape, search_shape,   &
       control, ndim_2d, "psmile_mg_method_2d_dble: eof")

#endif /* DEBUG_PRINT */
#endif /* DEBUG */

#ifdef VERBOSE
  print 9980, trim(ch_id), lev

  call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

9990 format (1x, a, ': psmile_mg_method_2d_dble: level', i3, &
             ', control', 6i6)
9980 format (1x, a, ': psmile_mg_method_2d_dble: eof, level', i3)

end subroutine psmile_mg_method_2d_dble