psmile_neigh_nearx_sub_irr_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_Neigh_nearx_sub_irr_dble
!
! !INTERFACE:

      subroutine psmile_neigh_nearx_sub_irr_dble ( grid_id,  &
                           x_coords, y_coords, z_coords,     &
                           coords_shape,                     &
                           mask_array, mask_shape, mask_available, &
                           sin_values, cos_values,           &
                           grid_valid_shape, search_mode,    &
                           neighbors_3d, num_neigh, nloc,    &
                           extra_search, ijk,                &
                           sin_search, cos_search, z_search, &
                           jbeg, jend, width, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_nearx_sub_irr_dble

      implicit none
!
! !DEFINED PARAMETERS:
!
!  NREF_3D = Number of locations in 3d to be controlled
!          = 4 ** ndim_3d
!  NREF_2D = Number of locations in 2d to be controlled
!          = 4 ** ndim_2d
!
!  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
!  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
!  dble_earth = Earth radius in meter
!
      Integer, Parameter              :: nref_3d = 4 * 4 * 4
      Integer, Parameter              :: nref_2d = 4 * 4
      Integer, Parameter              :: lon = 1
      Integer, Parameter              :: lat = 2
      Integer, Parameter              :: maxlen = 64
!
      Double Precision, Parameter     :: dble_earth = 6400000.0
!
      Double Precision, Parameter     :: tol = 1d-6 ! muss argument werden
      Double Precision, Parameter     :: eps = 1d-6
      Double Precision, Parameter     :: acosp1 =  1.0d0
      Double Precision, Parameter     :: acosm1 = -1.0d0
!
! !INPUT PARAMETERS:
!
      Integer,          Intent (In)   :: grid_id

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

      Integer,          Intent (In)   :: nloc
!
!     Total number of locations to be transferred
!
      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))
      Double Precision, Intent (In)   :: y_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2))
      Double Precision, Intent (In)   :: z_coords(coords_shape(1,3): 
                                                  coords_shape(2,3))

!     Coordinates of the method

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

!     Dimension of (method) mask array

      Logical,          Intent (In)   :: mask_array (mask_shape (1,1): 
                                                     mask_shape (2,1), 
                                                     mask_shape (1,2): 
                                                     mask_shape (2,2), 
                                                     mask_shape (1,3): 
                                                     mask_shape (2,3))
!     Mask of the method

      Logical,          Intent (In)   :: mask_available

!     Is mask specified in array "mask_array" ?
!     mask_available = .false. : Mask is not available
!                      .true.  : Mask is     available

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

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

      Double Precision, Intent (In)   :: sin_values (grid_valid_shape(1,1): 
                                                     grid_valid_shape(2,1), 
                                                     grid_valid_shape(1,2): 
                                                     grid_valid_shape(2,2), 
                                                     2)
!
!     Sin values of Longitudes and Latitudes (x_coords, y_coords)
!
      Double Precision, Intent (In)   :: cos_values (grid_valid_shape(1,1): 
                                                     grid_valid_shape(2,1), 
                                                     grid_valid_shape(1,2): 
                                                     grid_valid_shape(2,2), 
                                                     2)
!    
!     Cos values of Longitudes and Latitudes (x_coords, y_coords)
!
      Integer,          Intent (In)   :: search_mode
!
!     Specifies the search mode for nearest neigbours with
!        search_mode = 3 : Full 3d-search
!        search_mode = 2 : Search in 2d-hyperplane (1st 2 direction lon, lat)

      Integer,          Intent (In)   :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.
!
      Integer,          Intent (In)   :: jbeg, jend
!
!     Shape of input arrays
!
      Integer,          Intent (In)   :: ijk (ndim_3d, jbeg:jend)
!
!     Coarse grid index on level nlev-2 (for dimensions specified 
!     by search mode). Coarsening factor had to be 4.

      Double Precision, Intent (In)   :: sin_search (jbeg:jend, lat)
!
!     Sin values of the points for which the nearest neighbours
!     should be searched.
!
      Double Precision, Intent (In)   :: cos_search (jbeg:jend, lat)
!
!     Cos values of the points for which the nearest neighbours
!     should be searched.
!
      Double Precision, Intent (In)   :: z_search (jbeg:jend)
!
!     Z values of the points for which the nearest neighbours
!     should be searched.
!
      Integer,          Intent (In)   :: width (ndim_3d)
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Extra_search_info), Intent(InOut) :: extra_search
!
!     Number of locations where
!       (*) required mask values were not true
!
! !OUTPUT PARAMETERS:
!
      Integer,          Intent (Out)  :: neighbors_3d (ndim_3d, nloc, num_neigh)

!     Indices of neighbor locations (interpolation bases)
!
      Integer,          Intent (Out)  :: ierror

!     Returns the error code of PSMILE_Neigh_nearx_sub_irr_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
!     ... For locations controlled
!
      Integer                         :: nrefv (2:ndim_3d)
      Integer                         :: dim1 (2)
      Logical                         :: global_search
      Double Precision, Pointer       :: dist_dble (:, :)
!
!     ... (1d) Indices in (entire) data to be extra searched
!
      Integer                         :: i, j, k, jpart
      Integer                         :: len, leni, lenj, lenk, lenij
!
      Integer, Pointer                :: indices (:)
!
!     ... 3d Indices in local data to be extra searched
!
      Integer                         :: ii, jj, kk
      Integer                         :: n, nref, nrefd
!
      Integer                         :: ijkdst (ndim_3d, nref_3d)
      Integer                         :: ijkctl (ndim_3d, nref_3d)
      Integer                         :: irange (2, ndim_3d)
!
!     ... auxiliary arrays
!
      Integer                         :: itemp (ndim_3d)
!
!     ... for mask values
!
      Logical                         :: mask_dist (nref_3d)
      Logical                         :: mask_ind  (jbeg:jend)
!
!     ... for nearest neighbour control
!
      Integer                         :: nmin (1)
      Integer                         :: nfound
!
      Double Precision                :: dble_huge
      Double Precision                :: distance (nref_3d), dist, fac, val
      Double Precision                :: dist2 (4), dist3 (2*2*2)
      Double Precision                :: val2 (4)

      Integer                         :: nlev
      Type (Extra_search_dble)        :: arrays
      Type (Grid), Pointer            :: grid_info
      Type (Enddef_mg_double), Pointer :: my_arrays
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_nearx_sub_irr_dble" searches the next "num_neigh"
! nearest neighbours on the method-grid (x_coords, y_coords, z_coords)
! for the extra indices specified in "extra_search%indices".
!
! There are 2 different of search modes supported by the routine:
! search_mode = 3: The nearest neighbours are searched are searched
!                  within the 3d region.
! search_mode = 2: The nearest neighbours are searched only in the
!                  2d hyperplane (lon- and lat-direction).
!
! Subroutine "PSMILe_Neigh_nearx_sub_irr_dble" is designed for grids of
! of type "PRISM_Irrlonlat_regvrt".
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 17.02.05    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_nearx_sub_irr_dble.F90 2687 2010-10-28 15:15:52Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_nearx_sub_irr_dble.F90 2687 2010-10-28 15:15:52Z coquart $'
!
!----------------------------------------------------------------------
!
!  Relative control indices 
!      
! Schoener waere eine liste, die weniger sortierung nachher verlangt
! (dh. die nach den Summe der relativen Indices sortiert ist).
! Vielleicht sollte man immer annehmen, dass ijk in der "Mitte" liegt
!
      data ((ijkctl (i, n), i=1,ndim_3d), n = 1,nref_3d) &
            /  0, 0, 0,   1, 0, 0,   2, 0, 0,   3, 0, 0, &
               0, 1, 0,   1, 1, 0,   2, 1, 0,   3, 1, 0, &
               0, 2, 0,   1, 2, 0,   2, 2, 0,   3, 2, 0, &
               0, 3, 0,   1, 3, 0,   2, 3, 0,   3, 3, 0, &
               0, 0, 1,   1, 0, 1,   2, 0, 1,   3, 0, 1, &
               0, 1, 1,   1, 1, 1,   2, 1, 1,   3, 1, 1, &
               0, 2, 1,   1, 2, 1,   2, 2, 1,   3, 2, 1, &
               0, 3, 1,   1, 3, 1,   2, 3, 1,   3, 3, 1, &
               0, 0, 2,   1, 0, 2,   2, 0, 2,   3, 0, 2, &
               0, 1, 2,   1, 1, 2,   2, 1, 2,   3, 1, 2, &
               0, 2, 2,   1, 2, 2,   2, 2, 2,   3, 2, 2, &
               0, 3, 2,   1, 3, 2,   2, 3, 2,   3, 3, 2, &
               0, 0, 3,   1, 0, 3,   2, 0, 3,   3, 0, 3, &
               0, 1, 3,   1, 1, 3,   2, 1, 3,   3, 1, 3, &
               0, 2, 3,   1, 2, 3,   2, 2, 3,   3, 2, 3, &
               0, 3, 3,   1, 3, 3,   2, 3, 3,   3, 3, 3/
!
      data nrefv /nref_2d, nref_3d/
!                            
!----------------------------------------------------------------------
!
!  Initialization
!
!----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      grid_info => Grids(grid_id)
      indices   => extra_search%indices
      global_search = extra_search%global_search
!
      dble_huge = huge (distance(1))
      mask_ind  = .true.
!
!===> Allocate vector to store the distances found if required
!
      if (global_search) then
         dist_dble => extra_search%dist_dble
         dim1 (1) = 1
         dim1 (2) = extra_search%n_extra
      else
         dim1 (1) = jbeg
         dim1 (2) = jend
!
         Allocate (dist_dble (jbeg:jend, num_neigh), stat = ierror)

         if (ierror /= 0) then
            ierrp (1) = ierror
            ierrp (2) = num_neigh * (jend-jbeg+1)

            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'dist_dble', ierrp, 2, &
        __FILE__, __LINE__ )
            return
         endif
      endif
!
#ifdef PRISM_ASERTION
      if (global_search) then
         if (.not. Associated (extra_search%dist_dble)) then
            call psmile_assert (__FILE__, &
        __LINE__, 'extra_search%dist_dble was not allocated')
         endif
      endif
#endif

!
!===> Get number of reference points
!
      nrefd = nrefv (search_mode)
!
!===> Compute distances for each point in (jbeg:jend)
!
      do jpart = jbeg, jend, 5000

!cdir vector loopcnt=5000
      do j = jpart, min (jend, jpart+5000-1)
         i = indices(j)
!
!===> ... Get standard range to be controlled
!
         irange (1, :) = max (ijk(:, j),          grid_valid_shape (1, :))
         irange (2, :) = min (ijk(:, j)+width(:), 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 > nrefd) then
            print *, trim(ch_id), " nref, nredfd", nref, nrefd
            call psmile_assert (__FILE__, &
        __LINE__, 'nref > nrefd ')
         endif
#endif
!
!===> ... Set indices to be controlled
!
         if (nref == nrefd) then
               do n = 1, nrefd
               ijkdst (:, n) = ijk (:, j) + ijkctl (:, n)
               end do
         else
! 
            leni  =  irange (2, 1) - irange (1, 1) + 1
            lenij = (irange (2, 2) - irange (1, 2) + 1) * leni 
!
               do kk = 1, irange (2,3)-irange(1,3) + 1
               n = (kk-1)*lenij
                  do jj = 1, irange (2,2)-irange(1,2) + 1

!cdir vector
                     do ii = 1, leni
                     ijkdst (1, n+ii) = irange(1,1) + ii - 1
                     end do ! ii

                  ijkdst (2, n+1:n+leni) = irange (1,2) + jj - 1
                  n = n + leni
                  end do ! jj

               ijkdst (3, (kk-1)*lenij+1:kk*lenij) = irange (1,3) + kk - 1
               end do ! kk
         endif
!
!====> ... Get Mask values if necessary
!
         if (mask_available) then
!cdir vector
               do n = 1, nref
               mask_dist (n) = mask_array (ijkdst (1,n), ijkdst (2,n), &
                                           ijkdst (3,n))
               end do
!
            n = count (mask_dist(1:nref))
!
            if (n == 0) then
!              no value available
!
               nref = 0
               go to 10
            endif
!
            if (n < nref) then
               n = 0
                  do ii = 1, nref
                  if (mask_dist (ii)) then
                      n = n + 1
                      if (n /= ii) ijkdst (:, n) = ijkdst (:, ii)
                  endif
                  end do
               nref = n
            endif
         endif
!
!===> ... Compute distances to points
!
         fac = dble_earth + z_search (j)
!
!cdir vector
            do n = 1, nref
            val   = (  sin_values(ijkdst (1,n),                              &
                                  ijkdst (2,n), lat) * sin_search(j, lat)    &
                    +  cos_values(ijkdst (1,n),                              &
                                  ijkdst (2,n), lat) * cos_search(j, lat)    &
                    * (cos_values(ijkdst (1,n),                              &
                                  ijkdst (2,n), lon) * cos_search(j, lon)    &
                      +sin_values(ijkdst (1,n),                              &
                                  ijkdst (2,n), lon) * sin_search(j, lon)) )
#ifdef PRISM_ASSERTION
            if (val < acosm1 - eps .or. val > acosp1 + eps) then
               print *, i, j, val
               call psmile_assert (__FILE__, __LINE__, &
                                   'invalid value for acos')
            endif
#endif
            val = min (acosp1, val)
            val = max (acosm1, val)

            dist = fac * acos (val)

            distance (n) = dist*dist +                                 &
                           ( z_search (j)-z_coords(ijkdst (3,n)))**2
#ifdef DEBUGX
            if (i == 1387) then
               print *, 'dist, diff', dist, &
                        z_search (j)-z_coords(ijkdst (3,n))
            endif
#endif
            end do ! n
!
         nmin = minloc (distance(1:nref))
#ifdef MINLOCFIX
         if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
#ifdef DEBUGX
         if (i == 1387) then
            print *, 'psmile_neigh_nearx_sub_irr_dble: ictl =', i
            print *, 'ijkdst and distance:'
               do n = 1, nref
               print *, ijkdst (:, n), distance (n)
               end do
         endif
#endif
!
!===> ... Special case: Minimal distance == 0
!         Move other neighbours to "outside" of region.
!
         if (distance(nmin(1)) <= tol) then
            neighbors_3d (:, i, 1) = ijkdst (:, nmin(1))
!
               do n = 2, num_neigh
               neighbors_3d (:, i, n) = grid_valid_shape (1, :) - 1
               end do
!
            dist_dble (j, 1:num_neigh) = 0.0
!
            cycle
         endif
!
!===> ... Sort distances (silly sort)
!         Da muss noch etwas besseres her !
!
10       continue
!
         nfound = min (num_neigh, nref)
!
            do n = 1, nfound
!
            nmin = minloc (distance(n:nref))
#ifdef MINLOCFIX
            if (nmin(1).eq.0) nmin=1
#endif /* MINLOCFIX */
!
            if (nmin(1) /= 1) then
               jj = n + nmin (1) - 1
!
               itemp  (:)    = ijkdst (:, n)
               ijkdst (:, n) = ijkdst (:, jj)
               ijkdst (:, jj) = itemp  (:)
!
               dist = distance (n)
               distance (n) = distance (jj)
               distance (jj) = dist
            endif
            end do ! n
!
#ifdef PRISM_ASSERTION
            do n = 1, nfound-1
            if (distance(n) > distance (n+1)) exit
            end do
!
         if (n <= nfound-1) then
            print *, 'n =' , n
            print *, 'distance (n)  ', distance (n)
            print *, 'distance (n+1)', distance (n+1)
            print *, 'distance', distance (1:nfound)
            call psmile_assert (__FILE__, __LINE__, &
                                'incorrect sort')
         endif
#endif
!
!===> ... End of local search;
!         fill neighbors_3d and
!         store distances for global search if necessary
!
! Brauche ich die sqrt ?
!
         dist_dble (j, 1:nfound) = sqrt (distance (1:nfound))
         if (nfound < num_neigh) then
            dist_dble (j, nfound+1:num_neigh) = dble_huge

! das sollte man nicht machen (fuer nlev > 2)
! um mg suche einen startpunkt zu geben !
!
               do n = nfound+1, num_neigh
               ijkdst (:, n) = grid_valid_shape (1, :) - 1
               end do
         endif
!
         do n = 1, num_neigh
            neighbors_3d (:, i, n) = ijkdst (:, n)
         end do

! HR: Mein alter kommentar war falsch; man ist noch nicht fertig;
! Besser sollte man die lokale suche mit der originalen Zelle statt
! mit der groben Zelle machen.
! Trotzdem ist der jetzt gefundene Punkt schon eine gute
! Start-Loesung. Habe Derzeit keine Zeit weiter zu machen.
!        mask_ind(j) = .false.

      end do ! j
      end do ! jpart

      !
      !----------------------------------------------------------------------
      !
      ! 3rd) Do a wider search for those targets that are still missing data
      !
      !----------------------------------------------------------------------
      !
      if ( count(mask_ind) > 0 ) then
#define DEBUG_MG_EXTRA_SEARCH
#ifdef DEBUG_MG_EXTRA_SEARCH
         !-------------------------------------------------------------------
         !     Use a brut force for search in other regions
         !     A replacement for PSMILe_MG_srch_nneigh_irr_dble 
         !     as this requires further debugging
         !     Ref. ticket 47 with grid_medundemiT_cicle.nc as target
         !     and extra nneighbour search for bilinear interpolation. 
         !-------------------------------------------------------------------
         !
         call psmile_srch_nneigh_irreg2_dble (grid_id,   &
              search_mode, mask_ind,                     &
              sin_values, cos_values,                    &
              grid_valid_shape,                          &
              z_coords, coords_shape,                    &
              neighbors_3d, nloc, num_neigh,             &
              sin_search, cos_search, z_search,          &
              dist_dble, dim1, extra_search, jbeg, jend, &
              mask_array, mask_shape, mask_available,    &
              ierror)
#else
         !
         !----------------------------------------------------------------------
         !     Use multigrid info for search in other regions
         !
         !     What about points in boundary cells?
         !     For those we have to apply the additional search as well.
         !     corresponding to mask == .true. ?!!!!
         !----------------------------------------------------------------------
         !         
         if ( grid_info%nlev > 2 ) then
            !
            !===> Allocate memory and calculate distances on the coarsest level
            !     to begin the search for all points.
            !
            nlev = max(grid_info%nlev-3,2)
            !
            leni = grid_info%mg_infos(nlev)%levdim(1) + 1
            lenj = grid_info%mg_infos(nlev)%levdim(2) + 1
            lenk = grid_info%mg_infos(nlev)%levdim(3) + 1

            len  = leni*lenj*lenk

            ii = 0

            Allocate(arrays%radius(leni*lenj*lenk), STAT=jj); ii = ii + jj

            do kk = 1, 4
               Allocate(arrays%sin_ccorner_lon(kk)%vector(len), &
                        arrays%cos_ccorner_lon(kk)%vector(len), &
                        arrays%sin_ccorner_lat(kk)%vector(len), &
                        arrays%cos_ccorner_lat(kk)%vector(len), STAT=jj); ii = ii + jj
            enddo

            Allocate(arrays%sin_cmidp_lon%vector(len), &
                     arrays%cos_cmidp_lon%vector(len), &
                     arrays%sin_cmidp_lat%vector(len), &
                     arrays%cos_cmidp_lat%vector(len), STAT=jj); ii = ii + jj

            Allocate(arrays%sin_fmidp_lon%vector(maxlen), &
                     arrays%cos_fmidp_lon%vector(maxlen), &
                     arrays%sin_fmidp_lat%vector(maxlen), &
                     arrays%cos_fmidp_lat%vector(maxlen), STAT=jj); ii = ii + jj

            if ( ii /= 0 ) then
               ierrp (1) = ierror
               ierrp (2) = ii

               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'arrays%...%vector', ierrp, 2, &
                    __FILE__, __LINE__ )
               return

            endif

            my_arrays => grid_info%mg_infos(nlev)%double_arrays

            arrays%sin_ccorner_lon(1)%vector(1:len) = sin(my_arrays%chmin(1)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lon(1)%vector(1:len) = cos(my_arrays%chmin(1)%vector(1:len)*dble_deg2rad)
            arrays%sin_ccorner_lat(1)%vector(1:len) = sin(my_arrays%chmin(2)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lat(1)%vector(1:len) = cos(my_arrays%chmin(2)%vector(1:len)*dble_deg2rad)

            arrays%sin_ccorner_lon(2)%vector(1:len) = sin(my_arrays%chmin(1)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lon(2)%vector(1:len) = cos(my_arrays%chmin(1)%vector(1:len)*dble_deg2rad)
            arrays%sin_ccorner_lat(2)%vector(1:len) = sin(my_arrays%chmax(2)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lat(2)%vector(1:len) = cos(my_arrays%chmax(2)%vector(1:len)*dble_deg2rad)

            arrays%sin_ccorner_lon(3)%vector(1:len) = sin(my_arrays%chmax(1)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lon(3)%vector(1:len) = cos(my_arrays%chmax(1)%vector(1:len)*dble_deg2rad)
            arrays%sin_ccorner_lat(3)%vector(1:len) = sin(my_arrays%chmax(2)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lat(3)%vector(1:len) = cos(my_arrays%chmax(2)%vector(1:len)*dble_deg2rad)

            arrays%sin_ccorner_lon(4)%vector(1:len) = sin(my_arrays%chmax(1)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lon(4)%vector(1:len) = cos(my_arrays%chmax(1)%vector(1:len)*dble_deg2rad)
            arrays%sin_ccorner_lat(4)%vector(1:len) = sin(my_arrays%chmin(2)%vector(1:len)*dble_deg2rad)
            arrays%cos_ccorner_lat(4)%vector(1:len) = cos(my_arrays%chmin(2)%vector(1:len)*dble_deg2rad)

            arrays%sin_cmidp_lon%vector(1:len) = sin(my_arrays%midp(1)%vector(1:len)*dble_deg2rad)
            arrays%cos_cmidp_lon%vector(1:len) = cos(my_arrays%midp(1)%vector(1:len)*dble_deg2rad)
            arrays%sin_cmidp_lat%vector(1:len) = sin(my_arrays%midp(2)%vector(1:len)*dble_deg2rad)
            arrays%cos_cmidp_lat%vector(1:len) = cos(my_arrays%midp(2)%vector(1:len)*dble_deg2rad)

            n = 0

            if ( search_mode == 2 ) then
               !
               !===> ... 2d search mode
               !
               do k = 1, lenk
                  fac = dble_earth + my_arrays%midp(3)%vector(k)
                  !
                  do j = 1, lenj
                     !cdir vector
                     do i = 1, leni
                        n = n + 1

                        do ii = 1, 4
                           val2(ii) = ( arrays%sin_ccorner_lat(ii)%vector(n) *  &
                                        arrays%sin_cmidp_lat%vector(n)          &
                                      + arrays%cos_ccorner_lat(ii)%vector(n) *  &
                                        arrays%cos_cmidp_lat%vector(n)          &
                                      *(arrays%cos_ccorner_lon(ii)%vector(n) *  &
                                        arrays%cos_cmidp_lon%vector(n)          &
                                      + arrays%sin_ccorner_lon(ii)%vector(n) *  &
                                        arrays%sin_cmidp_lon%vector(n)) )
                           val2(ii) = min (acosp1, val2(ii))
                           val2(ii) = max (acosm1, val2(ii))

                        end do

                        dist2 = fac * acos (val2)
                        dist2 = dist2 * dist2

                        arrays%radius(n) = maxval (dist2)
                     enddo ! i
                  enddo
               enddo ! k

            else if ( search_mode == 3 ) then
               !
               !===> ... 3d search mode
               !
               do k = 1, lenk
                  fac = dble_earth + my_arrays%midp(3)%vector(k)
                  !
                  do j = 1, lenj
                     !cdir vector
                     do i = 1, leni
                        n = n + 1

                        do ii = 1, 4
                           val2(ii) = ( arrays%sin_ccorner_lat(ii)%vector(n) * &
                                        arrays%sin_cmidp_lat%vector(n)         &
                                      + arrays%cos_ccorner_lat(ii)%vector(n) * &
                                        arrays%cos_cmidp_lat%vector(n)         &
                                     * (arrays%cos_ccorner_lon(ii)%vector(n) * &
                                        arrays%cos_cmidp_lon%vector(n)         &
                                      + arrays%sin_ccorner_lon(ii)%vector(n) * &
                                        arrays%sin_cmidp_lon%vector(n)) )
                           val2(ii) = min (acosp1, val2(ii))
                           val2(ii) = max (acosm1, val2(ii))
                        end do

                        dist2 = fac * acos (val2)
                        dist2 = dist2 * dist2
                        !
                        dist3(1:4) = dist2 + (my_arrays%midp (3)%vector(k) - &
                                              my_arrays%chmin(3)%vector(k))**2

                        dist3(5:8) = dist2 + (my_arrays%midp (3)%vector(k) - &
                                              my_arrays%chmax(3)%vector(k))**2

                        arrays%radius(n) = maxval(dist3)

                     enddo ! i
                  enddo
               enddo ! k
            endif
            !
            ! werden die sqrt's wirklich gebraucht !?!?!
            ! Sehe ich nicht. a < b => sqrt (a) < sqrt (b)
            !
            arrays%radius(:) = sqrt (arrays%radius(:))
            !
            !===> ... Multigrid search
            !
            call psmile_mg_srch_nneigh_irr_dble (grid_id,   &
                 arrays, search_mode, nref_3d,              &
                 sin_values, cos_values, grid_valid_shape,  &
                 z_coords, coords_shape,                    &
                 neighbors_3d, nloc, num_neigh,             &
                 sin_search, cos_search, z_search,          & 
                 dist_dble, dim1, indices, jbeg, jend,      &
                 mask_ind,                                  &
                 mask_array, mask_shape, mask_available,    &
                 tol, ierror)
            if (ierror > 0) return
            !
            !===> ... Deallocate memory
            !
            Deallocate(arrays%sin_fmidp_lon%vector, arrays%cos_fmidp_lon%vector, &
                       arrays%sin_fmidp_lat%vector, arrays%cos_fmidp_lat%vector)

            Deallocate(arrays%sin_cmidp_lon%vector, arrays%cos_cmidp_lon%vector, &
                       arrays%sin_cmidp_lat%vector, arrays%cos_cmidp_lat%vector)

            do kk = 1, 4
               Deallocate(arrays%sin_ccorner_lon(kk)%vector, &
                          arrays%cos_ccorner_lon(kk)%vector, &
                          arrays%sin_ccorner_lat(kk)%vector, &
                          arrays%cos_ccorner_lat(kk)%vector)
            end do

            Deallocate(arrays%radius)

         endif ! grid_info%nlev > 2
#endif
      endif ! count(mask_ind) > 0

      if (.not. global_search) Deallocate (dist_dble)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      end subroutine PSMILe_Neigh_nearx_sub_irr_dble