psmile_mg_method_3d_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_method_3d_dble
!
! !INTERFACE:

      subroutine psmile_mg_method_3d_dble (                &
                           comp_info, nlev,                &
                           found, loc, range,              &
                           coords1, coords2, coords3,      &
                           shape, control,                 &
                           method_id,                      &
                           x_coords, y_coords, z_coords,   &
                           coords_shape,                   &
                           grid_valid_shape, cyclic,       &
                           chmin1, chmin2, chmin3,         &
                           chmax1, chmax2, chmax3,         &
                           midp1,  midp2,  midp3,          &
                           tol, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_mg_method_3d_dble

      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)   :: shape (2, ndim_3d)

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

      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_3d,               
                                                 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".

      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))
      Double Precision, Intent (In)   :: coords3 (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)   :: method_id

!     Method id

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

!     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), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))
      Double Precision, Intent (In)   :: y_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))
      Double Precision, Intent (In)   :: z_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))

!     Coordinates of the method

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

!     Specifies the valid block shape for the "ndim_3d"-dimensional block

!     Dimension of chmin, chmax and midp

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

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

      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),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      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),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      Double Precision, Intent (In)   :: chmin3 (grid_valid_shape(1,1)-1: 
                                                 grid_valid_shape(2,1),   
                                                 grid_valid_shape(1,2)-1: 
                                                 grid_valid_shape(2,2),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))

!     Minimum of grid coordinates per grid cell
!
! chmin (i, j, k) contains the minimal value of bounding box of
!              method cell ((i,j,k  ), (i+1,j,k  ), (i,j+1,k  ), (i+1,j+1,k  ),
!                           (i,j,k+1), (i+1,j,k+1), (i,j+1,k+1), (i+1,j+1,k+1))
!
! The first box
! chmin (range(1,1)-1,j,k) contains the minimal value of the bounding box
!                        of virtual cell
!       (array(range(1,1), j, k  ), array(range(1,1), j+1, k  ),
!        array(range(1,1), j, k+1), array(range(1,1), j+1, k+1),
!        bounding box of corner control volume of (range(1,1),j,k)
!
! The last box
! chmin (range(2,1),  j,k) contains the minimal value of the bounding box
!                        of virtual cell
!       (array(range(2,1), j, k  ), array(range(2,1), j+1, k  ),
!        array(range(2,1), j, k+1), array(range(2,1), j+1, k+1),
!        bounding box of corner control volume of (range(2,1),j,k)
!
! 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),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      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),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      Double Precision, Intent (In)   :: chmax3 (grid_valid_shape(1,1)-1: 
                                                 grid_valid_shape(2,1),   
                                                 grid_valid_shape(1,2)-1: 
                                                 grid_valid_shape(2,2),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))


!     Maximum of grid coordinates per grid cell

      Double Precision, Intent (In)   :: midp1  (grid_valid_shape(1,1)-1: 
                                                 grid_valid_shape(2,1),   
                                                 grid_valid_shape(1,2)-1: 
                                                 grid_valid_shape(2,2),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      Double Precision, Intent (In)   :: midp2  (grid_valid_shape(1,1)-1: 
                                                 grid_valid_shape(2,1),   
                                                 grid_valid_shape(1,2)-1: 
                                                 grid_valid_shape(2,2),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))
      Double Precision, Intent (In)   :: midp3  (grid_valid_shape(1,1)-1: 
                                                 grid_valid_shape(2,1),   
                                                 grid_valid_shape(1,2)-1: 
                                                 grid_valid_shape(2,2),   
                                                 grid_valid_shape(1,3)-1: 
                                                 grid_valid_shape(2,3))

!     Mid point of the cell

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

!     Returns the error code of PSMILE_mg_method_3d_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NREFD = Number of method coordinates to be controlled.
!        = 3 ** ndim_3d
!  NTET = Number of tetraeders used to determine whether
!         point is within the grid cell
!  EPS  = Tolerance for determinate of deformed grid cells
!  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.
!
#define USE6T
#ifdef USE6T
      Integer, Parameter              :: ntet = 6
#else
      Integer, Parameter              :: ntet = 4
#endif
      Double Precision, Parameter     :: eps = 1.0d-9
      Double Precision, Parameter     :: rmxeps = 1.0d30
      Double Precision, Parameter     :: dzero = 0.0
!
      Integer, Parameter              :: nrefd = 3 * 3 * 3
!
      Integer, Parameter              :: val_direct  =  1
      Integer, Parameter              :: val_coupler = -1
!
!     ... for levels
!
      Integer, Parameter              :: lev = 1
!
! !LOCAL VARIABLES
!
!     ... for tolerances
!
      Double Precision                :: tol1, tol2
!
!     ... for locations searched
!
      Integer                         :: i, j, k
      Integer                         :: ibegl, iendl, jpart, n
      Integer                         :: ic (ndim_3d)
!
!     ... For locations controlled
!
      Integer                         :: l, m
      Integer                         :: ijkctl (ndim_3d, nrefd)
      Integer                         :: ijkdst (ndim_3d, nrefd)
      Double Precision                :: dist (nrefd)
      Integer                         :: irange (2, ndim_3d)
!
      Integer                         :: ii, jj, kk, nref
      Integer                         :: nmin (1)
!   
!     ... Local reference point and edges
! 
      Integer                         :: ijkref (ndim_3d, ntet)
      Integer                         :: ijkedg (ndim_3d, 3, ntet)
!
      Double Precision                :: dble_huge
      Double Precision                :: xyzpnt (ndim_3d, ntet), 
                                         xyzref (ndim_3d, ntet)
      Double Precision                :: xyzedg (ndim_3d, 3, ntet)
      Double Precision                :: d1 (ntet), d2 (ntet), d3 (ntet)
      Double Precision                :: dnorm (ntet), det (ntet)
!
#ifdef DEBUG
!     ... for locations searched
!
      Integer                         :: nfnd (0:2)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_method_3d_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
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_method_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_method_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices
! ??? ist das richtig ??? muessen das nicht mehr sein
!
      data ((ijkctl (i, n), i=1,ndim_3d), n = 1,nrefd) &
                            / -1,-1,-1,   0,-1,-1,   1,-1,-1, &
                              -1, 0,-1,   0, 0,-1,   1, 0,-1, &
                              -1, 1,-1,   0, 1,-1,   1, 1,-1, &
                              -1,-1, 0,   0,-1, 0,   1,-1, 0, &
                              -1, 0, 0,   0, 0, 0,   1, 0, 0, &
                              -1, 1, 0,   0, 1, 0,   1, 1, 0, &
                              -1,-1, 1,   0,-1, 1,   1,-1, 1, &
                              -1, 0, 1,   0, 0, 1,   1, 0, 1, &
                              -1, 1, 1,   0, 1, 1,   1, 1, 1/
!
!  Relative indices of tetraheders
!
#ifdef USE6T
!         The cell is split into 6 tetraheders
!
      data ((ijkref (i, n), i = 1, ndim_3d), n = 1, ntet) &
           / 0, 0, 1,   0, 0, 1,   1, 0, 1,   0, 1, 1,    &
             0, 1, 1,   1, 1, 1 /
!
      data (((ijkedg (i,j, n), i = 1, ndim_3d), j = 1, 3), n = 1, ntet) &
           / 0, 0, 0,   1, 0, 0,   0, 1, 0,                             &
             1, 0, 0,   0, 1, 0,   1, 1, 0,                             &
             1, 0, 0,   1, 1, 0,   0, 0, 1,                             &
             0, 1, 0,   1, 1, 0,   0, 0, 1,                             &
             1, 1, 0,   0, 0, 1,   1, 0, 1,                             &
             1, 1, 0,   1, 0, 1,   0, 1, 1 /
#else
!         The cell is split into 4 tetraheders
!
      data  ((ijkref (i,n), i = 1, ndim_3d), n = 1, ntet)               &
           / 0, 0, 0,   1, 1, 0,   0, 1, 1,   1, 0, 1/
!
      data (((ijkedg (i,j,n), i = 1, ndim_3d), j = 1, 3), n = 1, ntet)  &
           / 1, 0, 0,   0, 1, 0,   0, 0, 1,                             &
             0, 1, 0,   1, 0, 0,   1, 1, 1,                             &
             1, 1, 1,   0, 0, 1,   0, 1, 0,                             &
             0, 0, 1,   1, 1, 1,   1, 0, 0 /
#endif /* USE6T */

!----------------------------------------------------------------------
!
!  Initialization
!
      ierror = 0

#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
!
! temporaeren speicherbereich (control (2,1)-control(1,1)+1)*ndtry anlegen
! und do n =1,nref vor i-loop ziehen ?!?
!
! ? statt immer alle nref zu kontrollieren, einmal location in der
! Zelle bestimmen ?
!
! ? ein-dimensionaler Loop statt i,j,k loop
!
      tol1 = tol + 1.0
      tol2 = tol * tol
!
      dble_huge = huge (dist(1))
!
#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)
!
!  Besser ein gatter/Scatter hier einbauen
!-----------------------------------------------------------------------
!
            ibegl = control (1, 1)
!
!===> ... Look for the next cells to be controlled
!
            do while (ibegl <= control(2,1))
!
!cdir vector
               do i = ibegl, control(2,1)
               if (found(i, j, k) == lev) exit
               end do
            if (i > control(2,1)) exit
            ibegl = i
!
!cdir vector
               do i = ibegl, control(2,1)
               if (found(i, j, k) /= lev) exit
               end do
!
            iendl = i - 1
!
!===> ... Look in the real cell (ic(1),ic(2),ic(3))
!
#ifdef PRISM_ASSERTION
               do i = ibegl, iendl
               if (loc(1, i,j,k) < coords_shape(1,1) .or. &
                   loc(1, i,j,k) > coords_shape(2,1) .or. &
                   loc(2, i,j,k) < coords_shape(1,2) .or. &
                   loc(2, i,j,k) > coords_shape(2,2) .or. &
                   loc(3, i,j,k) < coords_shape(1,3) .or. &
                   loc(3, i,j,k) > coords_shape(2,3)) exit
               end do
!
            if (i <= iendl) then
               print *, 'ic', loc (:, i,j,k)
               print *, 'coords_shape', coords_shape
               call psmile_assert (__FILE__, __LINE__, &
                                   'wrong index')
            endif
#endif
!
            do jpart = 0, (iendl-ibegl)/5000
!cdir vector, loopcnt=5000
loop_i:     do i = ibegl+jpart*5000, min (iendl, ibegl+(jpart+1)*5000-1)
!
            ic(:) = loc (:, i, j, k)
!
            irange (1, :) = max (ic(:)-1, grid_valid_shape (1, :))
            irange (2, :) = min (ic(:)+1, grid_valid_shape (2, :))
!
            nref = (irange (2,1) - irange (1,1) + 1) * &
                   (irange (2,2) - irange (1,2) + 1) * &
                   (irange (2,3) - irange (1,3) + 1)
!
#ifdef PRISM_ASSERTION
            if (nref <= 0) then
               call psmile_assert (__FILE__, __LINE__, &
                                   'nref <= 0')
            endif
#endif
!
            if (nref == nrefd) then
                  do n = 1, nrefd
                  ijkdst (:, n) = ic (:) + ijkctl (:, n)
                  end do
            else
                n = 0
                do kk = irange (1, 3), irange (2, 3)
                   do jj = irange (1, 2), irange (2, 2)
                      do ii = irange (1, 1), irange (2, 1)
                      n = n + 1
                      ijkdst (1, n) = ii
                      ijkdst (2, n) = jj
                      ijkdst (3, n) = kk
                      end do ! ii
                   end do ! jj
                end do ! kk
            endif
!
!===> ... Compute distance to points
!
               do n = 1, nref
               dist (n) = (x_coords (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                             coords1(i,j,k)) ** 2 + &
                          (y_coords (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                             coords2(i,j,k)) ** 2 + &
                          (z_coords (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                             coords3(i,j,k)) ** 2
               end do
!
!===> ... Look for the minimum distance to corner points of adjoining method cells
!
            nmin = MINLOC (dist(1:nref))
#ifdef MINLOCFIX
            if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
            if (dist(nmin(1)) < tol2) then
!
!===> ...... Point is located on the grid point of the method.
!
!              found (i,j,k)  = val_direct ! == lev == 1
               loc (:, i,j,k) = ijkdst (:, nmin(1))
!
#ifdef DEBUG
               nfnd (1) = nfnd (1) + 1
#endif /* DEBUG */

               cycle loop_i
            endif 
!
            found (i,j,k) = val_coupler
!
!===> ....... Get cell range "irange" for search
!
            if (ic(1) == grid_valid_shape(1,1) .or. &
                ic(2) == grid_valid_shape(1,2) .or. &
                ic(3) == grid_valid_shape(1,3)) then
               irange (1, :) = max (ic(:)-1, grid_valid_shape(1,:)-1)
               irange (2, :) = min (ic(:)+1, grid_valid_shape(2,:))
!
               nref = (irange (2,1) - irange (1,1) + 1) * &
                      (irange (2,2) - irange (1,2) + 1) * &
                      (irange (2,3) - irange (1,3) + 1)
!
#ifdef PRISM_ASSERTION
               if (nref <= 0) then
                  call psmile_assert (__FILE__, __LINE__, &
                                      'nref <= 0')
               endif
#endif
!
               if (nref == nrefd) then
                     do n = 1, nrefd
                     ijkdst (:, n) = ic (:) + ijkctl (:, n)
                     end do
               else
                   n = 0
                   do kk = irange (1, 3), irange (2, 3)
                      do jj = irange (1, 2), irange (2, 2)
                         do ii = irange (1, 1), irange (2, 1)
                         n = n + 1
                         ijkdst (1, n) = ii
                         ijkdst (2, n) = jj
                         ijkdst (3, n) = kk
                         end do ! ii
                      end do ! jj
                   end do ! kk
               endif
            endif
!
!===> ... Compute distance to the mid-points
!
!cdir vector
              do n = 1, nref
              if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n)) .and.   &
                  coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
                                           ijkdst(2,n),         &
                                           ijkdst(3,n))) then

                  dist (n) = (midp1 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                                coords1(i,j,k)) ** 2 + &
                             (midp2 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                                coords2(i,j,k)) ** 2 + &
                             (midp3 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                                coords3(i,j,k)) ** 2
               else
                  dist (n) = dble_huge
               endif
               end do
!
!-------------------------------------------------------------------------------
!     ...... Control whether the point is located within the cell starting in
!            "ic"
!            ??? Kann man das vereinfachen, wenn man weiss, dass der Punkt
!                in der bounding box ist ?
!-------------------------------------------------------------------------------
!
!===> ... Look for the minimum distance
!
            nmin = MINLOC (dist(1:nref))
#ifdef MINLOCFIX
            if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
                  do while (dist(nmin(1)) < dble_huge) 
                  ic = ijkdst (:, nmin(1))
!
                  if (ic(1) >= grid_valid_shape(1,1) .and. &
                      ic(1) <  grid_valid_shape(2,1) .and. &
                      ic(2) >= grid_valid_shape(1,2) .and. &
                      ic(2) <  grid_valid_shape(2,2) .and. &
                      ic(3) >= grid_valid_shape(1,3) .and. &
                      ic(3) <  grid_valid_shape(2,3)) then
!
!===> ....... Point is not located in a virtual method cell.
!             Get reference points of tetraeders for method cell (ic:ic+1)
!
                     do n = 1, ntet
                     xyzref (1, n) = x_coords (ic(1)+ijkref(1, n), &
                                               ic(2)+ijkref(2, n), &
                                               ic(3)+ijkref(3, n))
                     xyzref (2, n) = y_coords (ic(1)+ijkref(1, n), &
                                               ic(2)+ijkref(2, n), &
                                               ic(3)+ijkref(3, n))
                     xyzref (3, n) = z_coords (ic(1)+ijkref(1, n), &
                                               ic(2)+ijkref(2, n), &
                                               ic(3)+ijkref(3, n))
                     end do
!
                  xyzpnt (1, :) = coords1 (i, j, k) - xyzref (1, :)
                  xyzpnt (2, :) = coords2 (i, j, k) - xyzref (2, :)
                  xyzpnt (3, :) = coords3 (i, j, k) - xyzref (3, :)
!
!===> ....... Compute edges of tetraeders
!
                     do n = 1, ntet
                        do m = 1, 3
                        xyzedg (1, m, n) = x_coords (ic(1)+ijkedg(1, m, n),   &
                                                     ic(2)+ijkedg(2, m, n),   &
                                                     ic(3)+ijkedg(3, m, n)) - &
                                           xyzref (1, n)
                        xyzedg (2, m, n) = y_coords (ic(1)+ijkedg(1, m, n),   &
                                                     ic(2)+ijkedg(2, m, n),   &
                                                     ic(3)+ijkedg(3, m, n)) - &
                                           xyzref (2, n)
                        xyzedg (3, m, n) = z_coords (ic(1)+ijkedg(1, m, n),   &
                                                     ic(2)+ijkedg(2, m, n),   &
                                                     ic(3)+ijkedg(3, m, n)) - &
                                           xyzref (3, n)
                        end do
                     end do
!
!===> ........ Solve trilinear system; compute spat product
!
!  ?!?! Normierung, damit det nicht zu klein wird ?!
!
!   (xyzedg(*,1), xyzedg(*,2), xyzedg(*,3)) (x) = xyzpnt
!
!   det = a1*(b2*c3-b3*c2) + b1*(c2*a3-c3*a2) + c1*(a2*b3-a3*b2)
!
! :'a,'bs/c\([1-3]\)/xyzedg(\1,3,n)/g
! :'a,'bs/xyzedg(\([1-3]\),3,/xyzpnt(\1,  /g
!
!
                  det(:) = xyzedg(1,1,:)*(xyzedg(2,2,:)*xyzedg(3,3,:) - &
                                          xyzedg(3,2,:)*xyzedg(2,3,:))  &
                         + xyzedg(1,2,:)*(xyzedg(2,3,:)*xyzedg(3,1,:) - &
                                          xyzedg(3,3,:)*xyzedg(2,1,:))  &
                         + xyzedg(1,3,:)*(xyzedg(2,1,:)*xyzedg(3,2,:) - &
                                          xyzedg(3,1,:)*xyzedg(2,2,:))
!
                     do n = 1, ntet
                     dnorm (n) = dzero
                        do m = 1, 3
!cdir vector
                           do l = 1, ndim_3d
                           dnorm (n) = dnorm(n) + abs(xyzedg(l,m,n))
                           end do
                        end do
                     end do
!
#define ALLOW_DEGRADED_CELL
#ifdef ALLOW_DEGRADED_CELL
!
! ??? ist das koscher ?? wird d1(n), d2(n), d3(n) == 0
!cdir vector
                     do n = 1, ntet
                     if (abs(det(n)) <= eps*dnorm(n)*dnorm(n)*dnorm(n) &
                         .or. dnorm(n) == dzero) then
                        det (n) = sign (rmxeps, det(n))
                     else
                        det (n) = 1.0d0 / det (n)
                     endif
                     end do
#else 
!cdir vector
                     do n = 1, ntet
                     if (abs(det(n)) <= eps*dnorm(n)*dnorm(n)*dnorm(n) &
                         .or. dnorm(n) == dzero) go to 170
                     end do
!
                  det (:) = 1.0d0/det (:)

#endif /* ALLOW_DEGRADED_CELL */
!
                  d1(:) = (xyzpnt(1,  :)*(xyzedg(2,2,:)*xyzedg(3,3,:) -  &
                                          xyzedg(3,2,:)*xyzedg(2,3,:))   &
                         + xyzedg(1,2,:)*(xyzedg(2,3,:)*xyzpnt(3,  :) -  &
                                          xyzedg(3,3,:)*xyzpnt(2,  :))   &
                         + xyzedg(1,3,:)*(xyzpnt(2,  :)*xyzedg(3,2,:) -  &
                                          xyzpnt(3,  :)*xyzedg(2,2,:)))  &
                        * det (:)
!
                  d2(:) = (xyzedg(1,1,:)*(xyzpnt(2  ,:)*xyzedg(3,3,:) -  &
                                          xyzpnt(3,  :)*xyzedg(2,3,:))   &
                         + xyzpnt(1,  :)*(xyzedg(2,3,:)*xyzedg(3,1,:) -  &
                                          xyzedg(3,3,:)*xyzedg(2,1,:))   &
                         + xyzedg(1,3,:)*(xyzedg(2,1,:)*xyzpnt(3,  :) -  &
                                          xyzedg(3,1,:)*xyzpnt(2,  :)))  &
                        * det (:)
!
                  d3(:) = (xyzedg(1,1,:)*(xyzedg(2,2,:)*xyzpnt(3,  :) -  &
                                          xyzedg(3,2,:)*xyzpnt(2,  :))   &
                         + xyzedg(1,2,:)*(xyzpnt(2,  :)*xyzedg(3,1,:) -  &
                                          xyzpnt(3,  :)*xyzedg(2,1,:))   &
                         + xyzpnt(1,  :)*(xyzedg(2,1,:)*xyzedg(3,2,:) -  &
                                          xyzedg(3,1,:)*xyzedg(2,2,:)))  &
                        * det (:)
!
#ifdef USE6T
!
                     do n = 1, ntet 
                     if (min(d1(n), d2(n), d3(n)) >= -tol .and. &
                         max(d1(n),dzero) +                     &
                         max(d2(n),dzero) +                     &
                         max(d3(n),dzero) <=  tol1) go to 160
                     end do
!
                  go to 190
#else
                     do n = 1, ntet
                     if (min(d1(n), d2(n), d3(n)) < -tol .or. &
                         max(d1(n), d2(n), d3(n)) >  tol1) go to 190
                     end do
                  go to 160
#endif
                  endif
!
!===> ... point is contained in cell "ic"
!
160               continue
#ifdef DEBUG
                  nfnd (2) = nfnd (2) + 1
#endif /* DEBUG */
!
                  loc (:, i, j, k) = ic
!
                  cycle loop_i
!
!===> ... cell is degraded
!
#ifndef ALLOW_DEGRADED_CELL
170               print 9960, trim(ch_id), method_id, ic, &
                     x_coords(ic(1),   ic(2),   ic(3)),  &
                     y_coords(ic(1),   ic(2),   ic(3)),  &
                     z_coords(ic(1),   ic(2),   ic(3)),  &
                     x_coords(ic(1)+1, ic(2),   ic(3)),  &
                     y_coords(ic(1)+1, ic(2),   ic(3)),  &
                     z_coords(ic(1)+1, ic(2),   ic(3)),  &
                     x_coords(ic(1),   ic(2)+1, ic(3)),  &
                     y_coords(ic(1),   ic(2)+1, ic(3)),  &
                     z_coords(ic(1),   ic(2)+1, ic(3)),  &
                     x_coords(ic(1)+1, ic(2)+1, ic(3)),  &
                     y_coords(ic(1)+1, ic(2)+1, ic(3)),  &
                     z_coords(ic(1)+1, ic(2)+1, ic(3)),  &
                     x_coords(ic(1),   ic(2),   ic(3)+1), &
                     y_coords(ic(1),   ic(2),   ic(3)+1), &
                     z_coords(ic(1),   ic(2),   ic(3)+1), &
                     x_coords(ic(1)+1, ic(2),   ic(3)+1), &
                     y_coords(ic(1)+1, ic(2),   ic(3)+1), &
                     z_coords(ic(1)+1, ic(2),   ic(3)+1), &
                     x_coords(ic(1),   ic(2)+1, ic(3)+1), &
                     y_coords(ic(1),   ic(2)+1, ic(3)+1), &
                     z_coords(ic(1),   ic(2)+1, ic(3)+1), &
                     x_coords(ic(1)+1, ic(2)+1, ic(3)+1), &
                     y_coords(ic(1)+1, ic(2)+1, ic(3)+1), &
                     z_coords(ic(1)+1, ic(2)+1, ic(3)+1)
!
                     do n = 1, ntet
                     if (abs(det(n)) < eps) then
                        print 9950, n, det (n), dnorm(n), &
                          xyzedg (1, 1, n),                   &
                          xyzedg (1, 2, n), xyzedg (1, 3, n), &
                          xyzedg (2, 1, n), xyzedg (2, 2, n), &
                          xyzedg (2, 3, n), xyzedg (3, 1, n), &
                          xyzedg (3, 2, n), xyzedg (3, 3, n)
                     endif
!
                     end do
#endif /* ALLOW_DEGRADED_CELL */
!
!===> ... Point was not found in cell
!         Look for the minimum distance of others cells
!
190               continue
!
                  dist(nmin(1)) = dble_huge
                  nmin = MINLOC (dist(1:nref))
#ifdef MINLOCFIX
                  if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
               end do ! while (dist(nmin(1)) < dble_huge)
!
!===> ... Point was not found in any cell
!         Point may be located in a neighbouring block
!
!         Get nearest cell mid point (weiss nicht, ob das so clever ist ??)
!cdir vector
              do n = 1, nref
              dist (n) = (midp1 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                          coords1(i,j,k)) ** 2 + &
                         (midp2 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                          coords2(i,j,k)) ** 2 + &
                         (midp3 (ijkdst(1,n), ijkdst (2,n), ijkdst (3,n))-&
                          coords3(i,j,k)) ** 2
              end do

              nmin = MINLOC (dist(1:nref))
#ifdef MINLOCFIX
              if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
              loc (:, i,j,k) = ijkdst (:, nmin(1))
!
#ifdef DEBUG
               nfnd (0) = nfnd (0) + 1

#ifdef DEBUGX
               print *, trim(ch_id), ' Not found: -----'
               print *, 'i,j,k :', i, j, k
               print *, 'coords:', coords1(i,j,k), coords2(i,j,k), coords3(i,j,k)
               print *, 'irange:', irange
               print *, 'chmin low:', chmin1 (irange(1,1), irange(1,2), irange(1,3)), &
                                      chmin2 (irange(1,1), irange(1,2), irange(1,3)), &
                                      chmin3 (irange(1,1), irange(1,2), irange(1,3))
               print *, 'chmax low:', chmax1 (irange(1,1), irange(1,2), irange(1,3)), &
                                      chmax2 (irange(1,1), irange(1,2), irange(1,3)), &
                                      chmax3 (irange(1,1), irange(1,2), irange(1,3))
               print *, 'chmin hig:', chmin1 (irange(2,1), irange(2,2), irange(2,3)), &
                                      chmin2 (irange(2,1), irange(2,2), irange(2,3)), &
                                      chmin3 (irange(2,1), irange(2,2), irange(2,3))
               print *, 'chmax hig:', chmax1 (irange(2,1), irange(2,2), irange(2,3)), &
                                      chmax2 (irange(2,1), irange(2,2), irange(2,3)), &
                                      chmax3 (irange(2,1), irange(2,2), irange(2,3))
#endif /* DEBUGX */
    
#endif /* DEBUG */
!
            end do loop_i ! i
            end do ! jpart
!
!===> ... Is there an cyclic index
!
            do l = 1, ndim_3d
            if (cyclic(l)) then
!
!              This is an cyclic index;
!              Move points, which are located in virtual cell 0
!              into cell "grid_valid_shape(2,l)"
!
!cdir vector
                  do i = ibegl, iendl
                  if (loc (l, i,j,k) < grid_valid_shape (1,l)) then
                      loc (l, i,j,k) = grid_valid_shape (2,l)
                  endif
                  enddo
#ifdef REMOVED
            else
!
!              Move points located in virtual cell 0 into
!              cell 1; wird derzeit nur gemacht, um Probleme
!              in anderen Routinen zu verhindern
!              spaeter: diese indices in extra_search aufnehmen
!              oder nearest neighbour
!cdir vector
                  do i = ibegl, iendl
                  loc (l, i,j,k) = max (loc(l, i,j,k), grid_valid_shape (1,l))
                  enddo
#endif
!
            endif
            end do ! l
!
!===> ... Start index of the next search
!
            ibegl = iendl + 2
            end do ! while
!
            end do ! j
         end do ! k
!
!===> All done
!
#ifdef DEBUG
      print 9970, trim(ch_id), lev, nfnd
9970  format (1x, a, ': PSMILe_mg_method_3d_dble: lev =', i3, &
              ': not :', i5, ', dir: ', i6, ', fnd :', i6, i5)
#endif /* DEBUG */

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

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

9990 format (1x, a, ': psmile_mg_method_3d_dble: level', i3, &
                    ', control', 6i6)
9980 format (1x, a, ': psmile_mg_method_3d_dble: eof, level', i3)
9960 format (1x, a, ': psmile_mg_method_3d_dble: method id =', i5, &
            /1x,  ' Grid cell (', i4, ',', i4, ',', i4,           &
                 ') is degraded! Cell coordinates are:', 1p,      &
            8 (/1x, '(', 2 (e17.9, ','), e17.9, ')'))
9950 format (1x, i1, '-th spawn; det =', 1p, e17.9, ' (', e17.9, &
             ') of', 3 (/1x, 3e17.9, ','))


      end subroutine psmile_mg_method_3d_dble