psmile_mg_srch_nneigh_irr_real.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_srch_nneigh_irr_real
!
subroutine psmile_mg_srch_nneigh_irr_real (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_real, dim1, indices, jbeg, jend,      &
                mask_ind,                                  &
                mask_array, mask_shape, mask_available,    &
                tol, ierror)
!
! !USES:
!
  use PRISM_constants
!
  use PSMILe, dummy_interface => PSMILe_MG_srch_nneigh_irr_real

  implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id
!
!     Info on the component in which the donor cells
!     should be searched.
!
      Type (Extra_search_real)        :: arrays
!
!     Memory for metric information
!
      Integer, Intent (In)            :: nref_3d
!
!     Something like the maximum number of neighbor to be searched for
!
      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 2nd direction lon, lat)
!        search_mode = 1 : Search in 1d-plane
!
      Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
!
!     Specifies the valid block shape for the "ndim_3d"-dimensional block

      Real, 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)
!
      Real, 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)            :: coords_shape (2, ndim_3d)

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

      Real, Intent (In)               :: z_coords(coords_shape(1,3): 
                                                  coords_shape(2,3))

      Integer, Intent (In)            :: nloc
!
!     Second dimension of neighbors array "neighbors_3d" and
!     Total number of locations to be transferred
!
      Integer, Intent (In)            :: num_neigh
!
!     Number of neighbors to be searched.
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.
!
      Integer,          Intent (In)   :: jbeg, jend
!
!     Shape of input arrays
!
      Integer, Parameter              :: lat = 2
      Real, Intent (In)               :: sin_search (jbeg:jend, lat)
!
!     Sin values of the points for which the nearest neighbours
!     should be searched.
!
      Real, Intent (In)               :: cos_search (jbeg:jend, lat)
!
!     Cos values of the points for which the nearest neighbours
!     should be searched.
!
      Real, Intent (In)               :: z_search (jbeg:jend)
!
!     Z values of the points for which the nearest neighbours
!     should be searched.
!
      Integer, Intent (In)            :: dim1 (2)
!
!     Dimensions of "dist_real" in first direction 
!
! Rene: without interface description declaration as indices(:)
!       produces a segmentation fault. In this case it is
!       necessary to define indices from jbeg to jend.
!       This however violates assumptions for global search.
!
      Integer, Intent (In)            :: indices (:)
!
!     Indices of the extra points in entire list of all points to be
!     searched
!
      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
!
      Logical, Intent (In)             :: mask_ind(jbeg:jend)
!
!     Mask for j-indices already treated
!
      Real, Intent (In)                :: tol
!
!     Some tolerance may be used later to evaluate distances
!
! !INPUT/OUTPUT PARAMETERS:
!
      Real, Intent (InOut)             :: dist_real (dim1(1):dim1(2), num_neigh)
!
!     Initial distances and final distances
!
      Integer,           Intent (Out) :: neighbors_3d (ndim_3d, nloc, num_neigh)

!     Indices of neighbor locations (interpolation bases)
!     Index array for points that have been found
!
! !OUTPUT PARAMETERS:
!
     Integer, Intent (Out)           :: ierror

!    Returns the error code of PSMILe_MG_srch_nneigh_irr_real;
!            ierror = 0 : No error
!            ierror > 0 : Severe error!
!
! !DEFINED PARAMETERS:
!
!  real_earth = Earth radius in meter
!
  Real, Parameter                 :: real_earth = 6400000.0
!
  Real, Parameter                 :: acosp1 =  1.0
  Real, Parameter                 :: acosm1 = -1.0

  Integer, Parameter              :: lon = 1
!
! !LOCAL VARIABLES
!
  Real                            :: dist, dist1, distmax, val, fac, fac2

  Real                            :: sin_corner_lat(4,8)
  Real                            :: cos_corner_lat(4,8)
  Real                            :: sin_corner_lon(4,8)
  Real                            :: cos_corner_lon(4,8)

  Real                            :: sin_midp_lat(8)
  Real                            :: cos_midp_lat(8)
  Real                            :: sin_midp_lon(8)
  Real                            :: cos_midp_lon(8)

  Real                            :: distc (4)

  Real, Allocatable               :: distrad_c(:), radius_c(:)
  Real, Allocatable               :: distrad_f(:,:), radius_f(:,:)

  Logical, Allocatable            :: maskrad_c(:)
  Logical, Allocatable            :: maskrad_f(:,:)

  Integer, Allocatable            :: ijk0 (:,:)
  Integer                         :: ijk  (2,ndim_3d)
  Integer                         :: ii, jj, kk, nn
  Integer, Allocatable            :: next(:,:)
  Integer                         :: i, j, k, n
  Integer                         :: ifix, jfix, kfix
  Integer                         :: ncount
  Integer                         :: iloc(1)

  Type (Grid), Pointer            :: grid_info
  Type (Enddef_mg_real), Pointer :: my_arrays
!
!     ... multi-grid control
!
  Integer                         :: level, nlev
  Integer, Allocatable            :: boxsize (:, :)
!
!     ... number of points in i, j, k directions
!
  Integer                         :: len, leni, lenij, lenijk, lenj, lenk
  Integer                         :: leni_level, lenij_level
  Integer                         :: len_alloc
!
! jind = loop index in compact vector "sin_search", "cos_search",
!        "z_search" and "dist_real"
!
! ind  = Corresponding index in "neighbors_3d"
!
  Integer                         :: ind, jind
!
!     ... for error parameters
!
  Integer, parameter              :: nerrp = 1
  Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_srch_nneigh_irr_real" searches the next "num_neigh"
! nearest neighbours on the method-grid. num_neigh = 1 for the moment only.
!
! TODO: Ausnutzen der einkommenden Info "neighbors_3d" ueber die Indices
!       der schon gefundenen Punkte
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 05/02/01    R. Redler    created
!
!EOP
!----------------------------------------------------------------------
!
#ifdef VERBOSE
  print 9990, trim(ch_id)

  call psmile_flushstd
#endif /* VERBOSE */
!
!----------------------------------------------------------------------
!  Initialization
!----------------------------------------------------------------------
!
  ierror = 0

  grid_info => Grids(grid_id)

  nlev = grid_info%nlev - 3

  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

  lenij  = leni*lenj
  if ( search_mode == 3 ) then
     lenijk = lenij*lenk
  else
     lenijk = lenij
  endif
!
!----------------------------------------------------------------------
!  Allocate memory
!----------------------------------------------------------------------
!
  allocate (maskrad_c(lenijk), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = lenijk
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'maskrad_f', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  allocate (distrad_c(lenijk), radius_c(lenijk), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = lenijk
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'distrad_c, radius_c', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  if ( search_mode == 2 ) then
      len_alloc = 64
  else
      len_alloc = 64*8
  endif

  allocate (maskrad_f(len_alloc,nlev), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = 2*ndim_3d*nlev
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'maskrad_f', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  allocate (distrad_f(len_alloc,nlev), radius_f(64,nlev), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = 6*2*ndim_3d*nlev
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'distrad_f, radius_f', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  allocate (boxsize(ndim_3d,nlev), ijk0(ndim_3d,nlev), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = 2*ndim_3d*nlev
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'boxsize, ijk0', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  allocate (next(ndim_3d,nlev), stat = ierror)
  if (ierror /= 0) then
     ierrp (1) = ndim_3d*nlev
     ierror = PRISM_Error_Alloc
     call psmile_error ( ierror, 'next', ierrp, 1, &
          __FILE__, __LINE__ )
     return
  endif

  ijk0(:,nlev) = grid_info%ijk0
!
! compute coarsening factor to get from level "lev" onto the
! finest level "1". boxsize stores the number of data points
! in each direction in a particular box on a particular level.
!
  boxsize (1:ndim_3d, 1) = 1

  do level = 2, nlev
     do i = 1, ndim_3d
        if ( grid_info%mg_infos(level  )%levdim(i) /= &
             grid_info%mg_infos(level-1)%levdim(i) ) then
           boxsize(i, level) = boxsize(i, level-1) * 2
        else
           boxsize(i, level) = boxsize(i, level-1)
        end if
     end do
  end do

  do jind = jbeg, jend

     if ( mask_ind(jind) ) then
        !
        ! Get number of points in i,j,k direction on the start level
        !
        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

        lenij  = leni*lenj
        if ( search_mode == 3 ) then
           lenijk = lenij*lenk
        else
           lenijk = lenij
        endif

        maskrad_c = .true.

        ind = indices (jind)

        if ( search_mode == 2 ) then
           !
           ! Find nearest level in the vertical
           ! To do: Replace stupid search
           !
           dist = huge(dist)
           do k = grid_valid_shape(1,3), grid_valid_shape(2,3)
              val = abs( z_search (jind) - z_coords(k) )
              if ( val < dist ) then
                 val = dist
                 kfix = k
              endif
           enddo
        endif

        distmax = maxval (dist_real (jind, 1:num_neigh) )
        !
        ! distance from point to mid point of control volume on start level 
        ! Is this starting on level nlev too expensive (-> lenijk !)
        !
        fac = real_earth + z_search (jind)

        !cdir vector
        do n = 1, lenijk
           val = (  arrays%sin_cmidp_lat%vector(n) * sin_search(jind, lat)    &
                +   arrays%cos_cmidp_lat%vector(n) * cos_search(jind, lat)    &
                *  (arrays%cos_cmidp_lon%vector(n) * cos_search(jind, lon)    &
                +   arrays%sin_cmidp_lon%vector(n) * sin_search(jind, lon)) )

           val = min (acosp1, val)
           val = max (acosm1, val)

           dist = fac * acos (val)
           distrad_c (n) = dist * dist
        enddo

        if ( search_mode == 3 ) then
           my_arrays => grid_info%mg_infos(nlev)%real_arrays
           n = 0
           do k = 1, lenk
              distrad_c (n+1:n+lenij) = distrad_c (n+1:n+lenij) + &
                   (my_arrays%midp(3)%vector(k)-z_search(jind))**2
              n = n + lenij
           enddo
        endif

        ! distance from point to sphere
        ! Note: This distance may be negative if point is inside the volume.

        distrad_c (1:lenijk) = sqrt (distrad_c(1:lenijk)) - arrays%radius(1:lenijk)

200     continue ! entry point to proceed with other coarse level boxes

        ! reset pointers and values to fit coarse level.

        my_arrays => grid_info%mg_infos(nlev)%real_arrays

        next = 0

        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

        lenij  = leni*lenj
        if ( search_mode == 3 ) then
           lenijk = lenij*lenk
        else
           lenijk = lenij
        endif
        !
        ! ignore control volumes that are farther away than what we have already
        !
        do i = 1, lenijk
           maskrad_c (i) = maskrad_c(i) .and. (distrad_c(i) < distmax)
        enddo
        !
        ! if we are already done take next target point
        !
        if ( count(maskrad_c) < 1 ) cycle

        ncount = 0

        do i = 1, lenijk

           if ( count(maskrad_c) < 1 ) exit

           iloc(:) = minloc (distrad_c(1:lenijk), maskrad_c(1:lenijk))
#ifdef MINLOCFIX
           if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */

           kk = (iloc(1)                - 1) / lenij + 1
           jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
           ii =  iloc(1) - (kk-1)*lenij - (jj-1)*leni

           ijk(1,1) = ijk0(1,nlev) + (ii-1)*boxsize(1,nlev)
           ijk(1,2) = ijk0(2,nlev) + (jj-1)*boxsize(2,nlev)

           ijk(2,1) = min (grid_valid_shape(2,1), ijk(1,1)+boxsize(1,nlev)-1)
           ijk(2,2) = min (grid_valid_shape(2,2), ijk(1,2)+boxsize(2,nlev)-1)

           if ( search_mode == 2 ) then
              ijk(1,3) = kfix
              ijk(2,3) = kfix
           else
              ijk(1,3) = ijk0(3,nlev) + (kk-1)*boxsize(3,nlev)
              ijk(2,3) = min (grid_valid_shape(2,3), ijk(1,3)+boxsize(3,nlev)-1)
           endif

           ncount = count( mask_array(ijk(1,1):ijk(2,1), &
                                      ijk(1,2):ijk(2,2), &
                                      ijk(1,3):ijk(2,3)) )
           if ( ncount > 0 ) exit

           maskrad_c(iloc(1)) = .false. ! flag coarse cell as visited

        enddo

        if ( count(maskrad_c) < 1 ) cycle

        maskrad_c(iloc(1)) = .false.

        level = nlev

        maskrad_f = .true.

100     continue ! start f-cycle

        ! sphere around the control volume for all but the start level

        if ( level < nlev ) then
           !
           !----------------------------------------------------------------------
           !       Prepare next level
           !----------------------------------------------------------------------
           !
           ! Offset of local range ijk  relative to the level 1 index
           !
           leni   = boxsize(1,level+1)/boxsize(1,level)
           lenj   = boxsize(2,level+1)/boxsize(2,level)
           lenk   = boxsize(3,level+1)/boxsize(3,level)

           lenij  = leni*lenj
           if ( search_mode == 3 ) then
              lenijk = lenij*lenk
           else
              lenijk = lenij
           endif
           !
           ! Offset for "level" in multigrid array
           ! next(:level) origin is (0,0,0)
           !
           next(1,level) = (next(1,level+1)+ii-1) * leni
           next(2,level) = (next(2,level+1)+jj-1) * lenj
           next(3,level) = (next(3,level+1)+kk-1) * lenk
           !
           ! Corresponding offset for full domain
           !
           ijk0(1,level) = next(1,level)
           ijk0(2,level) = next(2,level)
           ijk0(3,level) = next(3,level)

           do i = level-1, 1, -1
              ijk0(1,level) = ijk0(1,level)*boxsize(1,i+1)/boxsize(1,i)
              ijk0(2,level) = ijk0(2,level)*boxsize(2,i+1)/boxsize(2,i)
              ijk0(3,level) = ijk0(3,level)*boxsize(3,i+1)/boxsize(3,i)
           enddo

           ijk0(1,level) = ijk0(1,level) + ijk0(1,nlev)
           ijk0(2,level) = ijk0(2,level) + ijk0(2,nlev)
           ijk0(3,level) = ijk0(3,level) + ijk0(3,nlev)
           !
           ! Point to next finer level
           !
           my_arrays => grid_info%mg_infos(level)%real_arrays

           leni_level  = grid_info%mg_infos(level)%levdim(1)+1
           lenij_level = leni_level * (grid_info%mg_infos(level)%levdim(2)+1)
           !
           ! restrict lengths to stay within valid bounds at outer boxes
           !
           if ( next(1,level)+leni > grid_info%mg_infos(level)%levdim(1)+1 ) &
                leni = grid_info%mg_infos(level)%levdim(1) + 1 - next(1,level)

           if ( next(2,level)+lenj > grid_info%mg_infos(level)%levdim(2)+1 ) &
                lenj = grid_info%mg_infos(level)%levdim(2) + 1 - next(2,level)

           if ( next(3,level)+lenk > grid_info%mg_infos(level)%levdim(3)+1 ) &
                lenk = grid_info%mg_infos(level)%levdim(3) + 1 - next(3,level)

           lenij  = leni*lenj
           if ( search_mode == 3 ) then
              lenijk = lenij*lenk
           else
              lenijk = lenij
           endif
           !
           ! collect the geographical data ...
           !
           do k = 1, lenk
              do j = 1, lenj
                 do i = 1, leni

                    n  = (k-1)*lenij + (j-1)*leni    + i

                    nn = (next(3,level) + k-1)*lenij_level + &
                         (next(2,level) + j-1)*leni_level  + &
                          next(1,level) + i

                    sin_corner_lon(1,n) = sin(my_arrays%chmin(lon)%vector(nn)*real_deg2rad)
                    cos_corner_lon(1,n) = cos(my_arrays%chmin(lon)%vector(nn)*real_deg2rad)
                    sin_corner_lat(1,n) = sin(my_arrays%chmin(lat)%vector(nn)*real_deg2rad)
                    cos_corner_lat(1,n) = cos(my_arrays%chmin(lat)%vector(nn)*real_deg2rad)

                    sin_corner_lon(2,n) = sin(my_arrays%chmax(lon)%vector(nn)*real_deg2rad)
                    cos_corner_lon(2,n) = cos(my_arrays%chmax(lon)%vector(nn)*real_deg2rad)
                    sin_corner_lat(2,n) = sin(my_arrays%chmin(lat)%vector(nn)*real_deg2rad)
                    cos_corner_lat(2,n) = cos(my_arrays%chmin(lat)%vector(nn)*real_deg2rad)

                    sin_corner_lon(3,n) = sin(my_arrays%chmax(lon)%vector(nn)*real_deg2rad)
                    cos_corner_lon(3,n) = cos(my_arrays%chmax(lon)%vector(nn)*real_deg2rad)
                    sin_corner_lat(3,n) = sin(my_arrays%chmax(lat)%vector(nn)*real_deg2rad)
                    cos_corner_lat(3,n) = cos(my_arrays%chmax(lat)%vector(nn)*real_deg2rad)

                    sin_corner_lon(4,n) = sin(my_arrays%chmin(lon)%vector(nn)*real_deg2rad)
                    cos_corner_lon(4,n) = cos(my_arrays%chmin(lon)%vector(nn)*real_deg2rad)
                    sin_corner_lat(4,n) = sin(my_arrays%chmax(lat)%vector(nn)*real_deg2rad)
                    cos_corner_lat(4,n) = cos(my_arrays%chmax(lat)%vector(nn)*real_deg2rad)

                    sin_midp_lon(n) = sin(my_arrays%midp(1)%vector(nn)*real_deg2rad)
                    cos_midp_lon(n) = cos(my_arrays%midp(1)%vector(nn)*real_deg2rad)
                    sin_midp_lat(n) = sin(my_arrays%midp(2)%vector(nn)*real_deg2rad)
                    cos_midp_lat(n) = cos(my_arrays%midp(2)%vector(nn)*real_deg2rad)

                 enddo
              enddo
           enddo
           !
           ! ... and calculate new distances
           !
           dist1 = 0.0

           do k = 1, lenk
              fac2 = real_earth + my_arrays%midp(3)%vector(k)

              if ( search_mode == 3 ) then
                 dist1 = max ( (my_arrays%midp (3)%vector(k) -     &
                                my_arrays%chmin(3)%vector(k))**2,  &
                               (my_arrays%midp (3)%vector(k) -     &
                                my_arrays%chmax(3)%vector(k))**2)
              endif
              !cdir vector
              do n = (k-1)*lenij+1, k*lenij

                 do nn = 1, 4
                    val = ( sin_corner_lat(nn,n) * sin_midp_lat(n)   &
                         +  cos_corner_lat(nn,n) * cos_midp_lat(n)   &
                         * (cos_corner_lon(nn,n) * cos_midp_lon(n)   &
                         +  sin_corner_lon(nn,n) * sin_midp_lon(n)) )

                    val = min (acosp1, val)
                    val = max (acosm1, val)

                    dist = fac2 * acos (val)

                    distc (nn) = sqrt(dist*dist + dist1)
                 enddo

                 radius_f(n,level) = maxval (distc)
              enddo

           enddo

           ! distance from point to mid point of control volume
           !  for all but the start level

           !cdir vector
           do n = 1, lenijk
              val = ( sin_midp_lat(n) * sin_search(jind,lat)   &
                   +  cos_midp_lat(n) * cos_search(jind,lat)   &
                   * (cos_midp_lon(n) * cos_search(jind,lon)   &
                   +  sin_midp_lon(n) * sin_search(jind,lon)) )

              val = min (acosp1, val)
              val = max (acosm1, val)

              dist = fac * acos (val)
              distrad_f(n,level) = dist * dist
           enddo

           if ( search_mode == 3 ) then
              n = 0
              do k = 1, lenk
                 distrad_f(n+1:n+lenij,level) = distrad_f (n+1:n+lenij,level) + &
                      (my_arrays%midp(3)%vector(k)-z_search(jind))**2
                 n = n + lenij
              end do
           endif

           distrad_f (1:lenijk,level) = sqrt (distrad_f(1:lenijk,level))

           ! distance from point to be searched for to sphere

           distrad_f(1:lenijk,level) = distrad_f(1:lenijk,level) - radius_f(1:lenijk,level)

           ! ignore control volumes that are farther away than what we have already

           do i = 1, lenijk
              maskrad_f(i,level) = maskrad_f(i,level) .and. (distrad_f(i,level) <= distmax)
           enddo

           ! go up one level and find another iloc there

           if ( count(maskrad_f(1:lenijk,level)) < 1 )  then

500            continue

               ! no more boxes to check, reset current mask and distance and go up

               maskrad_f(:,level) = .true.
               distrad_f(:,level) = huge(dist)

               level = level + 1

               if ( level+1 > nlev ) go to 200

               leni   = boxsize(1,level+1)/boxsize(1,level)
               lenj   = boxsize(2,level+1)/boxsize(2,level)
               lenk   = boxsize(3,level+1)/boxsize(3,level)

               lenij  = leni*lenj
               if ( search_mode == 3 ) then
                  lenijk = lenij*lenk
               else
                  lenijk = lenij
               endif

               if ( count(maskrad_f(1:lenijk,level)) < 1 ) go to 500

               iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
               if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
               kk = (iloc(1)                - 1) / lenij + 1
               jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
               ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

               go to 100

           else

              iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
              if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
           endif

           maskrad_f(iloc(1),level) = .false.

        endif ! (level < nlev)

300     continue
        !
        !-------------------------------------------------------------------------
        !    Ignore points masked out
        !-------------------------------------------------------------------------
        !
        if ( mask_available .and. iloc(1) > 0 ) then

           ! check whether we have valid points in this volume

           do i = 1, lenijk

              kk = (iloc(1)                - 1) / lenij + 1
              jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
              ii =  iloc(1) - (kk-1)*lenij - (jj-1)*leni

              ijk(1,1) = min (ijk0(1,level) + (ii-1)*boxsize(1,level), grid_valid_shape(2,1))
              ijk(1,2) = min (ijk0(2,level) + (jj-1)*boxsize(2,level), grid_valid_shape(2,2))

              ijk(2,1) = min (grid_valid_shape(2,1), ijk(1,1)+boxsize(1,level)-1)
              ijk(2,2) = min (grid_valid_shape(2,2), ijk(1,2)+boxsize(2,level)-1)

              if ( search_mode == 2 ) then
                 ijk(1,3) = kfix
                 ijk(2,3) = kfix
              else
                 ijk(1,3) = min (ijk0(3,level) + (kk-1)*boxsize(3,level), grid_valid_shape(2,3))
                 ijk(2,3) = min (grid_valid_shape(2,3), ijk(1,3)+boxsize(3,level)-1)
              endif

              ncount = count( mask_array(ijk(1,1):ijk(2,1), &
                                         ijk(1,2):ijk(2,2), &
                                         ijk(1,3):ijk(2,3)) )
              if ( ncount > 0 ) exit

              ! reject control volume without any valid points

              maskrad_f(iloc(1),level) = .false.

              if ( count(maskrad_f(1:lenijk,level)) > 0 ) then
                   iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
              else
                   iloc(1) = 0
                   exit
              endif

           enddo
           !
           ! ------------------------------------------------------------------------
           ! Special case if only one valid point is left over in the current volume.
           ! ------------------------------------------------------------------------
           !
           if ( ncount == 1 ) then

              do k = ijk(1,3),ijk(2,3)
                 do j = ijk(1,2),ijk(2,2)
                    do i = ijk(1,1),ijk(2,1)
                       if ( mask_array(i,j,k) ) then
                          ifix = i
                          jfix = j
                          kfix = k
                       endif
                    enddo
                 enddo
              enddo
              !
              !  Compute distance of point "neighbors_3d (:, ind, 1)"
              !
              val = ( sin_values(ifix, jfix, lat) * sin_search(jind, lat)    &
                   +  cos_values(ifix, jfix, lat) * cos_search(jind, lat)    &
                   * (cos_values(ifix, jfix, lon) * cos_search(jind, lon)    &
                   +  sin_values(ifix, jfix, lon) * sin_search(jind, lon)) )

              val = min (acosp1, val)
              val = max (acosm1, val)

              dist = fac * acos (val)

              if ( search_mode == 2 ) then
                 dist = dist*dist
              else
                 dist = dist*dist + ( z_search (jind)-z_coords(k) )**2
              endif

              !  Wrong for num_neigh > 1 !
              !  We need to keep old values which have
              !  been found already in this case.

              if ( sqrt ( dist ) < distmax ) then

                 neighbors_3d (1, ind, 1) = ifix
                 neighbors_3d (2, ind, 1) = jfix
                 neighbors_3d (3, ind, 1) = kfix

                 distmax =  sqrt ( dist )
                 dist_real(jind,1) = distmax

              endif

              maskrad_f(iloc(1),level) = .false.

!             iloc(1) = 0

700           continue

              maskrad_f(:,level) = .true.
              distrad_f(:,level) = huge(dist)

              level = level + 1

              if ( level+1 > nlev ) go to 200

              leni   = boxsize(1,level+1)/boxsize(1,level)
              lenj   = boxsize(2,level+1)/boxsize(2,level)
              lenk   = boxsize(3,level+1)/boxsize(3,level)

              lenij  = leni*lenj
              if ( search_mode == 3 ) then
                 lenijk = lenij*lenk
              else
                 lenijk = lenij
              endif

              do i = 1, lenijk
                 maskrad_f(i,level) = maskrad_f(i,level) .and. (distrad_f(i,level) <= distmax)
              enddo

              if ( count(maskrad_f(1:lenijk,level)) < 1 ) go to 700

              iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
              if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
              kk = (iloc(1)                - 1) / lenij + 1
              jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
              ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

              go to 300

           endif ! ncount == 1

        endif ! mask_available .and. iloc(1) > 0

        if ( iloc(1) > 0 ) then
            maskrad_f(iloc(1),level) = .false.
        else

800        continue

           maskrad_f(:,level) = .true.
           distrad_f(:,level) = huge(dist)

           level = level + 1

           if ( level+1 > nlev ) go to 200

           leni   = boxsize(1,level+1)/boxsize(1,level)
           lenj   = boxsize(2,level+1)/boxsize(2,level)
           lenk   = boxsize(3,level+1)/boxsize(3,level)

           lenij  = leni*lenj
           if ( search_mode == 3 ) then
              lenijk = lenij*lenk
           else
              lenijk = lenij
           endif

           do i = 1, lenijk
              maskrad_f(i,level) = maskrad_f(i,level) .and. (distrad_f(i,level) <= distmax)
           enddo

           if ( count(maskrad_f(1:lenijk,level)) < 1 ) go to 800

           iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
           if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
           kk = (iloc(1)                - 1) / lenij + 1
           jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
           ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

           go to 300

        endif

        level = level - 1

        if ( level == 3 ) go to 400 ! leave f-cycle and go to finest level 1

        go to 100 ! continue f-cycle

400     continue
        !
        ! ------------------------------------------------------------------
        ! go over all points of the final selected control volume on level 3
        ! and find the nearest points.
        ! ------------------------------------------------------------------
        !
        if ( iloc(1) > 0 ) then

           leni   = ijk(2,1) - ijk(1,1) + 1
           lenj   = ijk(2,2) - ijk(1,2) + 1
           lenk   = ijk(2,3) - ijk(1,3) + 1

           lenij  = leni*lenj
           lenijk = lenij*lenk

           ! Compute distance from point to potential neighbours

           do k = ijk(1,3), ijk(2,3)
              do j = ijk(1,2), ijk(2,2)
                 do i = ijk(1,1), ijk(2,1)

                    n  = (k-ijk(1,3))*lenij + (j-ijk(1,2))*leni + i-ijk(1,1)+1

                    val = ( sin_values(i, j, lat) * sin_search(jind, lat)    &
                         +  cos_values(i, j, lat) * cos_search(jind, lat)    &
                         * (cos_values(i, j, lon) * cos_search(jind, lon)    &
                         +  sin_values(i, j, lon) * sin_search(jind, lon)) )

                    val = min (acosp1, val)
                    val = max (acosm1, val)

                    dist = fac * acos (val)

                    distrad_f(n,level) = dist*dist

                 end do
              end do
           end do

           if ( search_mode == 3 ) then
              n = 0
              do k = 1, lenk
                 distrad_f(n+1:n+lenij,level) = distrad_f(n+1:n+lenij,level) + &
                       ( z_search (jind)-z_coords(k) )**2
                 n = n + lenij
              end do
           endif

           distrad_f (1:lenijk,level) = sqrt (distrad_f(1:lenijk,level))

           ! ignore invalid points

           if ( mask_available ) then

              n = 0
              do k = ijk(1,3), ijk(2,3)
                 do j = ijk(1,2), ijk(2,2)
                    do i = ijk(1,1), ijk(2,1)
                       n = n+1
                       maskrad_f(n,level) = mask_array(i,j,k)
                    enddo
                 enddo
              enddo

              do i = 1, lenijk
                 maskrad_f(i,level) = maskrad_f(i,level) .and. (distrad_f(i,level) <= distmax)
              enddo

              ncount = count(maskrad_f(1:lenijk,level))
              len = min(num_neigh,ncount)

           else

              do i = 1, lenijk
                 maskrad_f(i,level) = distrad_f(i,level) <= distmax
              enddo

              len = min(num_neigh, lenijk)

           endif
           !
           ! Find the remaining num_neigh nearest neighbours
           ! in this final level 3 control volume. 
           !
           do n = 1, len

              do nn = 1, lenijk

                 iloc(:) = minloc(distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
                 if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
                 kk = (iloc(1)                - 1) / lenij + 1
                 jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
                 ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

                 i = max(grid_valid_shape(1,1),ijk(1,1)+ii-1)
                 j = max(grid_valid_shape(1,2),ijk(1,2)+jj-1)

                 if ( search_mode == 2 ) then
                    k = kfix
                 else
                    k = max(grid_valid_shape(1,3),ijk(1,3)+kk-1)
                 endif

                 neighbors_3d (1, ind, n) = i
                 neighbors_3d (2, ind, n) = j
                 neighbors_3d (3, ind, n) = k

                 distmax = distrad_f(iloc(1),level)
                 dist_real (jind, n) = distmax

                 maskrad_f(iloc(1),level) = .false.
                 exit ! lopp nn and go for next neighbour

              enddo ! nn

           enddo ! n

           ! reset current mask and distance and go up to level 4

600        continue

           maskrad_f(:,level) = .true.
           distrad_f(:,level) = huge(dist)

           level = level + 1

           if ( level+1 > nlev ) go to 200

           leni   = boxsize(1,level+1)/boxsize(1,level)
           lenj   = boxsize(2,level+1)/boxsize(2,level)
           lenk   = boxsize(3,level+1)/boxsize(3,level)

           lenij  = leni*lenj
           if ( search_mode == 3 ) then
              lenijk = lenij*lenk
           else
              lenijk = lenij
           endif

           do i = 1, lenijk
              maskrad_f(i,level) = maskrad_f(i,level) .and. (distrad_f(i,level) <= distmax)
           enddo

           if ( count(maskrad_f(1:lenijk,level)) < 1 ) go to 600

           iloc(:) = minloc (distrad_f(1:lenijk,level), maskrad_f(1:lenijk,level))
#ifdef MINLOCFIX
           if (iloc(1).eq.0) iloc(1)=1
#endif /* MINLOCFIX */
           kk = (iloc(1)                - 1) / lenij + 1
           jj = (iloc(1) - (kk-1)*lenij - 1) / leni  + 1
           ii =  iloc(1) - (kk-1)*lenij - (jj-1) * leni

           go to 300

        endif

     endif ! mask_ind

  enddo ! jind
  !
  !===> All done
  !
  if (mask_available) then
     Deallocate (boxsize)
  endif

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

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

end subroutine PSMILe_MG_srch_nneigh_irr_real