psmile_mg_final_2d_dble.F90

Go to the documentation of this file.

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

      subroutine psmile_mg_final_2d_dble (grid_id, nlev,             &
                      chmin1, chmin2, chmax1, chmax2,                &
                      midp1,  midp2,  levdim,                        &
                      found, loc, range,                             &
                      coords1, coords2, shape, control,              &
                      corners1, corners2, corner_shape, nbr_corners, & 
                      tol, ierror)

!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_mg_final_2d_dble
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:

      Integer, Intent (In)            :: grid_id

!     Grid id
!
      Integer, Intent (In)            :: nlev

!     Total number of levels

      Integer, Intent (In)                        :: levdim (ndim_2d)

!     Dimensions of chmin, chmax and midp on level 1

      Double Precision, Intent (In)               :: chmin1 (0:levdim(1), 
                                                             0:levdim(2))
      Double Precision, Intent (In)               :: chmin2 (0:levdim(1), 
                                                             0:levdim(2))

!     Minimum of grid coordinates per grid cell on level 1

      Double Precision, Intent (In)               :: chmax1 (0:levdim(1), 
                                                             0:levdim(2))
      Double Precision, Intent (In)               :: chmax2 (0:levdim(1), 
                                                             0:levdim(2))

!     Maximum of grid coordinates per grid cell on level 1

      Double Precision, Intent (In)   ::  midp1 (0:levdim(1), 
                                                 0:levdim(2))
      Double Precision, Intent (In)   ::  midp2 (0:levdim(1), 
                                                 0:levdim(2))
!     Mid point of the cell

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

!     Dimension of loc and found

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

!     Dimension of coordinate arrays "coord1" and "coord2"
!     which contain the coordinates to be searched.
!     Indices of the grid cell in which the point was found.
!     The indices are relative to corner shape "corner_shape" !
!     I.e. the indices are already transformed into corner/method dimensions.

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

!     Coordinates to be searched

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

!     Input  value: index range in "coords" to be searched
!     Output value: index range in "coords" for which a cell was found

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

!     Dimension of corner coordinate arrays

      Integer, Intent (In)            :: nbr_corners

!     Number of corners

      Double Precision,   Intent (In) ::                 
         corners1 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

      Double Precision,   Intent (In) ::                 
         corners2 ( corner_shape(1,1):corner_shape(2,1), 
                    corner_shape(1,2):corner_shape(2,2), nbr_corners)

!     Corner coordinates

      Double Precision, Intent (In)   :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0

      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) = ??:
!     Otherwise, not contained in the corner 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))

!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_mg_final_2d_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NREFD = Number of corner coordinates to be controlled.
!        = 3 ** ndim_2d
!  NTET  = Number of triangles used to determine whether
!          point is within the grid cell
!  NDTRY = Number of neighbouring (fine) grid cells
!          to be controlled.
!  EPS   = Tolerance for determinate of deformed grid cells
!
      Integer, Parameter              :: ntet = 3
!
      Double Precision, Parameter     :: eps = 1.0d-9
      Double Precision, Parameter     :: rmxeps = 1.0d30
      Double Precision, Parameter     :: dzero = 0.0
!
      Integer, Parameter              :: nrefd = 3 * 3
      Integer, Parameter              :: ndtry = 3 * 3 - 1
!
      Integer, Parameter              :: temp_len = 256
!
!     ... for levels
!
      Integer, Parameter              :: lev = 1
      Integer, Parameter              :: not_found = -2
!
! !LOCAL VARIABLES
!
!     ... for tolerances
!
      Double Precision                :: tol1, tol2

!     ... for Grid

      Type (Grid), Pointer            :: grid_info
      Integer, Pointer                :: ijk0 (:)
!
!     ... for locations searched
!
      Integer                         :: i, j, k
      Integer                         :: ibegl, ibeglx, iendl, n
      Integer                         :: ic (temp_len, ndim_2d)
!
      Integer                         :: ii
!
!     ... For locations controlled
!
      Integer                         :: nmin (temp_len, 1)
      Integer                         :: ijkctl (ndim_2d, ndtry)
      Integer                         :: ijkdst (temp_len, ndim_2d, ndtry)
      Double Precision                :: dist (temp_len, ndtry)
!   
!     ... Local reference point and edges
! 
      Double Precision                :: dble_huge
!   
!     ... 
! 
      Integer                         :: nc, nc_part, nprev
      Integer                         :: list (temp_len*2)
      Integer                         :: part_list (temp_len)
      Logical                         :: fnd (temp_len)
!
!     ... for pole cells
!
      Logical                         :: control_pole_cell
      Integer                         :: n_pole
      Integer                         :: corner_dim (1)
      Integer                         :: ind (temp_len)
      Integer, Pointer                :: pole_array (:)
!
!     ... for search on coarser levels
!
      Integer                         :: loc2   (ndim_2d, temp_len)
      Integer                         :: newijk (ndim_2d, temp_len)
      Double Precision                :: xyz    (ndim_2d, temp_len)
!
#ifdef DEBUG
!     ... for locations searched
!
      Integer                         :: nfnd (0:4)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_final_2d_dble" searches the donor cells
! on the real (final) grid for the subgrid coords by the sending process.
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 13.03.07    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_final_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_final_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices for neighbour cells (excluding center cell)
!
      data ((ijkctl (i, n), i=1,ndim_2d), n = 1,ndtry) &
                  /-1,-1,   0,-1,   1,-1, &
                   -1, 0,           1, 0, &
                   -1, 1,   0, 1,   1, 1/
!
!----------------------------------------------------------------------
!
!  Initialization
!
      ierror = 0

      grid_info => Grids(grid_id)
      ijk0 => grid_info%ijk0

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

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (tol < 0.0) then
         call psmile_assert (__FILE__, __LINE__, &
                             "tol must be >= 0.0")
      endif
#endif
!
      tol1 = tol + 1.0
      tol2 = tol * tol
!
      dble_huge = huge (tol1)
      iendl = control (2,1) + 1
!
!===> Have pole cells to be controlled ?
!
      control_pole_cell = grid_info%pole_covered .and. &
                          Associated(grid_info%corner_pointer%pole_array)
      if (control_pole_cell) then
         corner_dim (1) = corner_shape(2,1) - corner_shape(1,1) + 1
         pole_array => grid_info%corner_pointer%pole_array
         n_pole = Size (pole_array)
      endif
!
#ifdef DEBUG
      nfnd (:) = 0
#endif /* DEBUG */
!
#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
         i = ictl_ind(1)
         j = ictl_ind(2)
         k = control(1,3)

         if (found (i,j,k) == lev) then
            print 8990, ictl_ind (1:ndim_2d), &
                        found (i,j,k), loc (:, i,j,k)
         endif
8990  format (1x, '### final ictl_ind', 2i5, '; found, loc ', i3, 1x, 2i6)
      endif
#endif
!
!===>
!
         do k = control(1,3), control(2,3)
loop_j:     do j = control(1,2), control (2,2)
!
!-----------------------------------------------------------------------
!           Control location of point (coords1(i,j,k), coords2(i,j,k))
!           in real corner cell "loc(:, i,j,k)"
! list ueber i index ausdehnen ?
!-----------------------------------------------------------------------
!
! temporaeren speicherbereich (control (2,1)-control(1,1)+1)*ndtry anlegen
!
! ? ein-dimensionaler Loop statt i,j loop
!
            nprev = 0

            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))
!
               nc = 0
               i = ibegl
!
                  do while (nc < temp_len .and. i <= control(2,1))
                     ibegl = i
!cdir vector
                        do i = ibegl, min (ibegl+(temp_len*2-1)-nc, control(2,1))
                        if (found(i, j, k) == lev) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                        end do
                  end do
!
               if (nc == 0) exit loop_i ! exit while loop (ibegl < ...)
               if (nc < temp_len) then
                  ibegl = control(2,1) + 1
               else  if (nc > temp_len) then
                  nc = temp_len
                  ibegl = list (temp_len+1)
               else
                  ibegl = list (temp_len) + 1
               endif
!
!===> ... Look in the real cell (ic(i, 1),ic(i, 2))
!
            ic (1:nc, 1) = loc (1, list(1:nc), j, k)
            ic (1:nc, 2) = loc (2, list(1:nc), j, k)
!
!===> ...... Control whether the point "(list(1:nc),j,k)" is contained
!            in the corner cell "ic(1:nc,*)".
!
            call psmile_control_cell_2d_dble (grid_id,               &
                      ic, nc, temp_len, list, j, k,                  &
                      coords1,  coords2,  shape,                     &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, fnd, ierror)
            if (ierror > 0) return

#ifdef DEBUG_TRACE
            if (j == ictl_ind(2)) then
!              Is the point to be traced in the list of points ?
                  do i = 1, nc
                  if (list(i) == ictl_ind(1)) exit
                  end do
!
               if (i <= nc) then
                  print *, 'ictl: after 1st control cell'
                  print *, 'i,j, fnd, loc', list(i), j, fnd(i), &
                           loc (1, list(i), j, k), loc (2, list(i), j, k)
                  print *, 'coords', coords1(list(i), j, k), &
                                     coords2(list(i), j, k)
                  print *, 'chmin1', chmin1(loc (1, list(i), j, k)-ijk0(1),  &
                                            loc (2, list(i), j, k)-ijk0(2)), &
                                     chmax1(loc (1, list(i), j, k)-ijk0(1),  &
                                            loc (2, list(i), j, k)-ijk0(2))
                  print *, 'chmin2', chmin2(loc (1, list(i), j, k)-ijk0(1),  &
                                            loc (2, list(i), j, k)-ijk0(2)), &
                                     chmax2(loc (1, list(i), j, k)-ijk0(1),  &
                                            loc (2, list(i), j, k)-ijk0(2))

                     do n = 1, 4
                     print *, corners1 (loc (1, list(i), j, k),     &
                                        loc (2, list(i), j, k), n), &
                              corners2 (loc (1, list(i), j, k),     &
                                        loc (2, list(i), j, k), n)
                     enddo
               endif
            endif
#endif
!
!===> ... Store data on points which are contained in cell "ic"
!         Note: "loc" is not changed.
!
            n = count (fnd(1:nc))
#ifdef DEBUG
            nfnd (1) = nfnd (1) + n
#endif /* DEBUG */
!
            if (n == nc) cycle loop_i
!
!===> ... Get points which were not found
!
            nc_part = 0
!cdir vector
               do i = 1, nc
               if (.not. fnd (i)) then
                  nc_part =  nc_part + 1
                  part_list(nc_part) = list (i)
               endif
               end do
!
!===> ... Get compact list of points which were not found
!
            nc = nc_part
            list (1:nc) = part_list (1:nc)
!
!-----------------------------------------------------------------------
!           Control pole cells if necessary
!           Note: pole_array(n) contains the index of corner cell
!                 in 1-dimensional vector
!-----------------------------------------------------------------------
!
#define CONTROL_POLE_CELL
#ifdef CONTROL_POLE_CELL
            if (control_pole_cell) then
!cdir vector
                  do i = 1, nc
                  ind (i) = &
                    (loc (2,list(i),j,k)-corner_shape(1,2)) * corner_dim (1) + &
                     loc (1,list(i),j,k)-corner_shape(1,1) + 1
                  end do
!
               fnd (1:nc) = .false.
                  do n = 1, n_pole
                  fnd (1:nc) = fnd (1:nc) .or. pole_array(n) == ind (1:nc)
                  end do
!
               n = Count (fnd(1:nc))
#ifdef DEBUG
               nfnd (1) = nfnd (1) + n
#endif /* DEBUG */
!
               if (n == nc) cycle loop_i
               if (n > 0) then
!
!===> ... Get points which were not found
!
                  nc_part = 0
!cdir vector
                     do i = 1, nc
                     if (.not. fnd (i)) then
                        nc_part =  nc_part + 1
                        part_list(nc_part) = list (i)
                     endif
                     end do
!
!===> ... Get compact list of points which were not found
!
                  nc = nc_part
                  list (1:nc) = part_list (1:nc)
               endif
            endif
#endif
!
!-----------------------------------------------------------------------
!          Look in the neigbourhood of the fine corner cell
!          loc(:,list(i),j,k) for the points which were not found
!          in this corner cell.
!          Note: "ijkdst" is MG location in (0:levdim(:))
!-----------------------------------------------------------------------
!
                do n = 1, ndtry
                    do ii = 1, ndim_2d
                    ijkdst (1:nc, ii, n) = max (0,               &
                       min (loc(ii,list(1:nc),j,k)-ijk0(ii)+ijkctl(ii,n), &
                            levdim(ii)))
                    end do
!
!cdir vector
                    do i = 1, nc
                    if (coords1(list(i),j,k) >= chmin1(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords2(list(i),j,k) >= chmin2(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords1(list(i),j,k) <= chmax1(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n)) .and.   &
                        coords2(list(i),j,k) <= chmax2(ijkdst(i, 1,n),         &
                                                       ijkdst(i, 2,n))) then

                       dist (i, n) = (coords1(list(i),j,k) -                    &
                                      midp1(ijkdst(i, 1,n), ijkdst(i, 2,n)))**2 &
                                   + (coords2(list(i),j,k) -                    &
                                      midp2(ijkdst(i, 1,n), ijkdst(i, 2,n)))**2

                    else

                       dist (i, n) = dble_huge

                    endif
!
#ifdef DEBUG_TRACE
                    if (list(i) == ictl_ind(1) .and. j == ictl_ind(2)) then
                       print 8980, list(i), j,  n, &
                                   (ijkdst (i, ii, n), ii = 1, ndim_2d)
                       if (dist(i,n) == dble_huge) then
                          print *, 'outside of box'
                       else
                          print 8970, dist (i, n)
                       endif
8980 format (1x, 'ictl: neighbour', 2i6, ', n', i4, ', ijkdst', 2i6)
8970 format (1x, 'dist', 1p, d23.16)
                       print *, chmin1(ijkdst(i, 1,n), ijkdst(i, 2,n)), &
                                coords1(list(i),j,k), &
                                chmax1(ijkdst(i, 1,n), ijkdst(i, 2,n))
                       print *, chmin2(ijkdst(i, 1,n), ijkdst(i, 2,n)), &
                                coords2(list(i),j,k), &
                                chmax2(ijkdst(i, 1,n), ijkdst(i, 2,n))
                    endif
#endif
                    end do ! i
                 end do ! n
!
!===> ... Look for the minimum distance
!
50            continue
#if 0
              do i = 1, nc
              nmin(i,:) = MINLOC (dist(i, 1:ndtry))
              end do
!
#ifdef MINLOCFIX
              do i = 1, nc
              if (nmin(i, 1) == 0) nmin (i, 1) = 1
              end do

#endif /* MINLOCFIX */
#else
               nmin (1:nc, 1) = 1
!
                  do n = 2, ndtry
!cdir vector
                     do i = 1, nc
                     if (dist(i, n) < dist(i, nmin(i,1))) then
                        nmin (i,1) = n
                     endif
                     end do
                  end do
#endif
!
!     ... Mark points which were not found in a fine cell 
!         by negative number for further search
!
!cdir vector
                  do i = 1, nc
                  if (dist(i, nmin(i, 1)) == dble_huge) then
                     found (list(i),j,k) = - lev
                  endif
                  end do
!
!     ... Look for points which were found in a fine grid cell (box)
!
            nc_part = 0
!cdir vector
               do i = 1, nc
               if (dist(i, nmin(i, 1)) /= dble_huge) then
                  nc_part = nc_part + 1
                  part_list (nc_part) = list (i)
!
                  ic (nc_part, :) = ijk0 (1:ndim_2d) + ijkdst (i, :, nmin(i, 1))
               endif
               end do
!
            nprev = nprev + (nc - nc_part)
            if (nc_part == 0) cycle loop_i
!
!           Create compact list of remaining points
!           HR ? Use/Create an index instead of moving dist and ijkdst ?
!
               do n = 1, nc
               if (dist(n, nmin(n, 1)) == dble_huge) exit
               end do
!
            nc_part = n - 1
!
               do i = n+1, nc
               if (dist(i, nmin(i, 1)) /= dble_huge) then
                  nc_part = nc_part + 1
!
                  nmin (nc_part, 1) = nmin (i, 1)
                  dist (nc_part, :) = dist(i, :)
                  ijkdst (nc_part, :, :) = ijkdst (i, :, :)
               endif
               end do
!
#ifdef DEBUG_TRACE
            if (j == ictl_ind(2)) then
               do i = 1, nc_part
               if (part_list(i) == ictl_ind(1)) exit
               end do
!
               if (i <= nc_part) then
                  print *, 'ictl fine', part_list(i), j, ', nmin', nmin(i,1), &
                           'ic', ic(i, :)
               endif
            endif
#endif
!
            nc = nc_part
            call psmile_control_cell_2d_dble (grid_id,               &
                      ic, nc, temp_len, part_list, j, k,             &
                      coords1,  coords2,  shape,                     &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, fnd, ierror)
            if (ierror > 0) return

#ifdef DEBUG_TRACE
            if (j == ictl_ind(2)) then
!              Is the point to be traced in the list of points ?
                  do i = 1, nc
                  if (list(i) == ictl_ind(1)) exit
                  end do
!
               if (i <= nc) then
                  print *, '--- after 2. control cell'
                  print *, 'i,j, fnd, loc', part_list(i), j, fnd(i), &
                           loc (1, part_list(i), j, k), &
                           loc (2, part_list(i), j, k)
                  print *, 'coords', coords1(part_list(i), j, k), &
                                     coords2(part_list(i), j, k)
                  print *, 'ic', ic (i, :)
               end if
            endif
#endif
!
            nc_part = count (fnd(1:nc))
!cdir vector
               do i = 1, nc
               if (fnd(i)) then
!                 Point was contained in real cell "ic(i, :)"
!                 found  (part_list(i),j,k) = lev
                  loc (:, part_list(i),j,k) = ic (i, :)
               endif
               end do
!
#ifdef DEBUG
            nfnd (2) = nfnd (2) + nc_part
#endif /* DEBUG */
            if (nc_part == nc) cycle loop_i
!
!           Try to look in other neighbouring cells
!
            nc_part = 0
!cdir vector
               do i = 1, nc
               if (.not. fnd(i)) then
                  dist(i, nmin(i,1)) = dble_huge

                  nc_part = nc_part + 1
                  list (nc_part) = part_list (i)
               endif
               end do
!
!cdir vector
               do n = 1, nc
               if (fnd(n)) exit
               end do
!
            nc_part = n - 1
               do i = n+1, nc
               if (.not. fnd(i)) then
                  nc_part = nc_part + 1
                  dist   (nc_part, :)    = dist (i, :)
                  ijkdst (nc_part, :, :) = ijkdst (i, :, :)
               endif
               end do
!
            nc = nc_part
            go to 50
!
         end do loop_i ! while
!
!===> ... Have points to be controlled on a coarser level ?
!
         if (nprev == 0) cycle loop_j
!
!-----------------------------------------------------------------------
!          Some points were not found.
!          Try to interpolate the fine point from neighbouring points
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!     Some points were still not found.
!     Look again on the next coarser level whether the points
!     are located in other cells.
!-----------------------------------------------------------------------
!
         if (lev+1 < nlev) then
            ibeglx = control (1, 1)
!
!===> ... Look for the next points to be controlled
!         (i.e. look for points, which were found on simplified level 1 == lev).
!         Store "i" indices of points in list "list".
!
loop_3:        do while (ibeglx < iendl)
!
               nc = 0
               i = ibeglx
!
                  do while (nc < temp_len .and. i < iendl)
                     ibeglx = i
!cdir vector
                        do i = ibeglx, min (ibeglx+((temp_len*2-1)-nc), &
                                            control(2,1))
! das ist falsch
                        if (found(i, j, k) == -lev) then 
                           nc = nc + 1
                           list (nc) = i
                        endif
                        end do
!
                     if (nc < temp_len .and. &
                         i > control(2,1) .and. i < iendl) then
                        i = control(2,1)
! das ist falsch
                        if (found (i,j,k) == -lev) then
                           nc = nc + 1
                           list (nc) = i
                        endif
                     endif 
                  end do
!
               if (nc == 0) exit loop_3 ! exit while loop (ibeglx < ...)
!
               if (nc < temp_len) then
                  ibeglx = control(2,1) + 1
               else if (nc > temp_len) then
                  nc = temp_len
                  ibeglx = list (temp_len+1)
               else
                  ibeglx = list (temp_len) + 1
               endif
!
!===> ... Look in the coarser parts of cell loc(:,list(i),j,k)
!     !!! TODO: collect xyz and reduce number of subroutine calls
! 3d loop mit nprev ?
!
!cdir vector
                  do i = 1, nc
                  xyz (1, i) = coords1 (list(i),j,k)
                  xyz (2, i) = coords2 (list(i),j,k)
!
                  loc2 (:, i) = loc (:, list(i),j,k)-ijk0(1:ndim_2d)
                  end do
!
                  call psmile_mg_final_prev_2d_dble (grid_id, nlev,  &
                      loc2, xyz, part_list, newijk, nc,              &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, ierror)
                  if (ierror > 0) return
!cdir vector
                  do i = 1, nc
                  if (part_list(i) == 0) then
!
!                    Point was not found
!
                     found (list(i),j,k) = not_found
!
#ifdef DEBUGX
print *, 'not found A: ', list(i),j,k
print *, 'chmin ', chmin1 (list(i)-ijk0(1),j-ijk0(2)), chmin2 (list(i)-ijk0(1),j-ijk0(2))
print *, 'chmax ', chmax1 (list(i)-ijk0(1),j-ijk0(2)), chmax2 (list(i)-ijk0(1),j-ijk0(2))
print *, 'loc', loc (1, list(i),j,k), loc (2, list(i),j,k)
print *, 'coords', coords1 (list(i),j,k), coords2 (list(i),j,k)
print *, 'chmin ', chmin1 (loc (1, list(i),j,k)-ijk0(1), loc (2, list(i),j,k)-ijk0(2)), &
                   chmin2 (loc (1, list(i),j,k)-ijk0(1), loc (2, list(i),j,k)-ijk0(2))
print *, 'chmax ', chmax1 (loc (1, list(i),j,k)-ijk0(1), loc (2, list(i),j,k)-ijk0(2)), &
                   chmax2 (loc (1, list(i),j,k)-ijk0(1), loc (2, list(i),j,k)-ijk0(2))
print *, 'loc2', loc2 (1, nc), loc2 (2, nc)
print *, 'loc(i-1, j,k)', found (list(i)-1,j,k), loc(:, list(i)-1, j,k) 
print *, 'loc(i+1, j,k)', found (list(i)+1,j,k), loc(:, list(i)+1, j,k) 
#endif /* DEBUGX */

#ifdef DEBUG
                     nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
                  else
                     found  (list(i),j,k) = lev
                     loc (:, list(i),j,k) = ijk0(1:ndim_2d) + newijk (:, i)
#ifdef DEBUG
                     nfnd (4) = nfnd (4) + 1
#endif /* DEBUG */
                  endif
                  end do
!
               end do loop_3 ! while
            else
!              The coarsest level contains only 1 cell.
!              All points which are not found are outside.

!cdir vector
               do i = control (1,1), control (2,1)
               if (found (i,j,k) == -lev) then
                  found (i,j,k) = not_found
#ifdef DEBUG
                  nfnd (0) = nfnd (0) + 1
#endif /* DEBUG */
               endif
               end do
            endif
!
         end do loop_j ! 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
         i = ictl_ind(1)
         j = ictl_ind(2)
         k = control(1,3)

         if (found (i,j,k) == lev .or. found(i,j,k) == not_found) then
            print 8990, ictl_ind (1:ndim_2d), &
                        found (i,j,k), loc (:, i,j,k)
!
            print *, 'chmin ', &
                     chmin1 (loc(1,i,j,k)-ijk0(1), loc(2,i,j,k)-ijk0(2)), &
                     chmin2 (loc(1,i,j,k)-ijk0(1), loc(2,i,j,k)-ijk0(2))
            print *, 'coords', coords1 (i,j,k), coords2 (i,j,k)
            print *, 'chmax ', &
                     chmax1 (loc(1,i,j,k)-ijk0(1), loc(2,i,j,k)-ijk0(2)), &
                     chmax2 (loc(1,i,j,k)-ijk0(1), loc(2,i,j,k)-ijk0(2))
         endif
      endif
#endif
!
#ifdef DEBUG
      print 9970, trim(ch_id), lev, nfnd
#endif /* DEBUG */

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

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

9990  format (1x, a, ': psmile_mg_final_2d_dble: level', i3, &
                     ', control', 6i6)
9980  format (1x, a, ': psmile_mg_final_2d_dble: eof, level', i3)
#ifdef DEBUG
9970  format (1x, a, ': psmile_mg_final_2d_dble: lev =', i3, &
              ': not fnd', i5, ', fnd =', i6, 3i5)
#endif

      end subroutine psmile_mg_final_2d_dble