psmile_neigh_near_3d_irr3_dble.F90

Go to the documentation of this file.

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

      subroutine psmile_neigh_near_3d_irr3_dble (          &
                           grid_id,                        &
                           coords1, coords2, coords3,      &
                           x_coords, y_coords, z_coords,   &
                           coords_shape,                   &
                           mask_array, mask_shape,  mask_available,    &
                           sin_values_lon, cos_values_lon, &
                           sin_values_lat, cos_values_lat, &
                           grid_valid_shape, search_mode,  &
                           srcloc, nsearch,                &
                           nloc, nprev,                    &
                           neighbors_3d, num_neigh,        &
                           extra_search, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_near_3d_irr3_dble

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id

!     Link to the 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)            :: nprev
!
!     Number of locations for which neighbours were already searched.
!
      Integer, Intent (In)            :: nsearch
!
!     Number of locations for which neighbours should be searched in this
!     sweep.
!
      Integer, Intent (In)            :: srcloc (ndim_3d, nloc)

!     Indices of the grid cell in which the point was found.
!
      Double Precision, Intent (In)   :: coords1 (nloc)
      Double Precision, Intent (In)   :: coords2 (nloc)
      Double Precision, Intent (In)   :: coords3 (nloc)

!     Coordinates of the target points (which were searched and found)

      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))
      Double Precision, Intent (In)   :: y_coords(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_lon (grid_valid_shape(1,1): 
                                                         grid_valid_shape(2,1))
      Double Precision, Intent (In)   :: sin_values_lat (grid_valid_shape(1,2): 
                                                         grid_valid_shape(2,2))
!
!     Sin values of Longitudes and Latitudes (x_coords, y_coords)
!
      Double Precision, Intent (In)   :: cos_values_lon (grid_valid_shape(1,1): 
                                                         grid_valid_shape(2,1))
      Double Precision, Intent (In)   :: cos_values_lat (grid_valid_shape(1,2): 
                                                         grid_valid_shape(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)
!        search_mode = 1 : Search in 1d-plane
!
      Integer, Intent (In)            :: num_neigh
!
!     Last dimension of neighbors array "neighbors_3d" and
!     number of neighbors to be searched.
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Extra_search_info)        :: 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_near_3d_irr3_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !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
!  NREF_1D = Number of locations in 1d to be controlled
!          = 4 ** ndim_1d
!
!  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
!  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
!  REAL_EARTH = Earth radius
!
      Integer, Parameter              :: nref_3d = 4 * 4 * 4
      Integer, Parameter              :: nref_2d = 4 * 4
      Integer, Parameter              :: nref_1d = 4
      Integer, Parameter              :: lon = 1
      Integer, Parameter              :: lat = 2
!
      Double Precision, Parameter     :: dble_earth = 6400.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
!
! !LOCAL VARIABLES
!
!     ... For locations controlled
!
      Double Precision, Pointer       :: sin_search (:, :)
      Double Precision, Pointer       :: cos_search (:, :)
!
      Integer                         :: nrefv (1:ndim_3d)
!
!     ... (1d) Indices in (entire) data to be extra searched
!
      Integer                         :: ibeg, iend
      Integer                         :: leni, lenij
!
!     ... 3d Indices in local data to be extra searched
!
      Integer                         :: i, ipart
      Integer                         :: ii, jj, kk
      Integer                         :: n, nref, nrefd, sd
!
      Integer, Pointer                :: ijk (:, :)
      Integer                         :: ijk0 (ndim_3d)
      Integer                         :: ijkdst (ndim_3d, nref_3d)
      Integer                         :: ijkctl (ndim_3d, nref_3d)
      Integer                         :: irange (2, ndim_3d)
      Integer                         :: width (ndim_3d)
!
!     ... auxiliary arrays
!
      Integer                         :: itemp (ndim_3d)
!
!     ... for mask values
!
      Logical                         :: mask_dist (nref_3d)
!
!     ... for nearest neighbour control
!
      Integer                         :: nfound
      Integer                         :: nmin (1)
!
      Double Precision                :: dble_huge
      Double Precision                :: dist, fac, val
      Double Precision                :: distance (nref_3d)
#ifdef TODO
      Double Precision                :: dist_control
#endif
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_near_3d_irr3_dble" searches the next "num_neigh"
! nearest neighbours on the method-grid (x_coords, y_coords, z_coords)
! for the subgrid coords sent by the destination process.
!
! There are 3 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).
! search_mode = 1: The nearest neighbours are searched only in the
!                  1d plane
!
! Subroutine "PSMILe_Neigh_near_3d_irr3_dble" is designed for
! (*) grids of type "PRISM_Reglonlatvrt" and
! (*) source locations "srcloc" which are stored as
!     3-dimensional vector
!     containing the 3d indices of the cells found for the target points
!     "coords1, coords2, coords3" to be searched.

!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_near_3d_irr3_dble.F90 2868 2011-01-10 15:45:22Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_near_3d_irr3_dble.F90 2868 2011-01-10 15:45:22Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices 
!      
#ifdef OLD
      Integer, Parameter              :: nrefd = 3 * 3 * 3
!
      data ((ijkctl (i, n), i=1,ndim_3d), n = 1,nref_3d) &
                            /  0, 0, 0,                       &
                               1, 0, 0,   0, 1, 0,   0, 0, 1, &
                              -1, 0, 0,   0,-1, 0,   0, 0,-1, &
                              -1,-1, 0,   1,-1, 0,            &
                              -1, 1, 0,   1, 1, 0,            &
                              -1, 0,-1,   1, 0,-1,            &
                              -1, 0, 1,   1, 0, 1,            &
                               0,-1,-1,   0, 1,-1,            &
                               0,-1, 1,   0, 1, 1,            &
                              -1,-1,-1,   1,-1,-1,            &
                              -1, 1,-1,   1, 1,-1, &
                              -1,-1, 1,   1,-1, 1, &
                              -1, 1, 1,   1, 1, 1/
#endif
! schoener waere eine liste, die weniger sortierung (s.o.) nachher verlangt
! veielleicht 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_1d, nref_2d, nref_3d/
!                            
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      ibeg = nprev + 1
      iend = nprev + nsearch
!
#ifdef PRISM_ASSERTION
!
!===> Internal control
!
      if (search_mode < 1 .or. search_mode > ndim_3d) then
         print *, trim(ch_id), " search_mode ", search_mode
         call psmile_assert (__FILE__, __LINE__, &
                             'search_mode is out of range ')
      endif
!
      if (nloc < iend ) then
         print *, trim(ch_id), " nloc, nprev, nsearch ", nloc, nprev, nsearch
         call psmile_assert (__FILE__, __LINE__, &
                             'nloc < nprev + PRODUCT (nlocs) ')
      endif
!
!===> Are the locations within the correct shape ?
!
!cdir vector
      do i = ibeg, iend
!
         if (srcloc(1,i) < coords_shape(1,1) - 1 .or. &
             srcloc(1,i) > coords_shape(2,1)     .or. &
             srcloc(2,i) < coords_shape(1,2) - 1 .or. &
             srcloc(2,i) > coords_shape(2,2)     .or. &
             srcloc(3,i) < coords_shape(1,3)     .or. &
             srcloc(3,i) > coords_shape(2,3)) exit
      end do
!
      if (i <= iend) then
         print *, "i, ibeg, iend", i, ibeg, iend
         print *, "srcloc ", srcloc(:, i)
         print *, "shape  ", coords_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index in srcloc')
      endif
#endif
!
      dble_huge = huge (distance(1))
!
!===> Get number of reference points
!
      nrefd = nrefv (search_mode)
!
!-----------------------------------------------------------------------
!     Allocate and set temporary array which is used to transform
!     the coordinates of poinst to be searched.
!     (*) Longitudes and Latitudes are transformed
!         from degrees into radients.
!     (*) The z-values are currently not transformed
!         ??? Whats about PRISM_Irrlonlat_sigmavrt, ...
!
!     Is it possible to save the values, in order to be reused.
!-----------------------------------------------------------------------
!
      Allocate (sin_search(ibeg:iend, lat), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = nsearch * lat

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'sin_search', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      Allocate (cos_search(ibeg:iend, lat), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = nsearch * lat

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'cos_search', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif

      Allocate (ijk(ndim_3d, ibeg:iend), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = nsearch * ndim_3d

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'ijk', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!===> Compute values required to compute distances
!
!cdir vector
      do i = ibeg, iend
         sin_search (i, lon) = coords1 (i) * dble_deg2rad
      end do
!
!cdir vector
      do i = ibeg, iend
         sin_search (i, lat) = coords2 (i) * dble_deg2rad
      end do
!
      cos_search = cos (sin_search)
      sin_search = sin (sin_search)
!
!-----------------------------------------------------------------------
!     Compute distances
!
!     dist = (dble_earth + coords3)*ACOS ( SIN(lat_coords)*SIN(lat_search) +
!                                          COS(lat_coords)*COS(lat_search) *
!                                         (COS(lon_coords)*COS(lon_search) +
!                                          SIN(lon_coords)*SIN(lon_search)) )
!     dist = SQRT(dist**2 + (coords3-z_coords) ** 2)
!
!     with xxx_search = Coordinates to be searched
!     and  xxx_coords = Grid (Methods) coordinates (transformed)
!
! Suchen mit 
!                    nlev = Grids(grid_id)%nlev
!                    Grids(grid_id)%mg_infos(lev)%double_arrays%chmin, &
!                    Grids(grid_id)%mg_infos(lev)%double_arrays%chmax, &
!
!-----------------------------------------------------------------------
!
      ijk0 = Grids(grid_id)%ijk0
      sd = search_mode
!
      width (1:sd) = 3
!
!===> Get coarse grid index on level nlev-2
!
      do i = 1, sd
         ijk (i, ibeg:iend) = ijk0(i) + &
                              ((srcloc(i, ibeg:iend) - ijk0(i)) / 4) * 4
      end do ! i
!
      do i = sd+1, ndim_3d
            ijk (i, ibeg:iend) = srcloc (i, ibeg:iend)
            width (i) = 0
      end do ! i
!
!===> Compute distances for each point in (ibeg:iend)
!
      do ipart = ibeg, iend, 5000
!cdir vector loopcnt=5000
      do i = ipart, min(iend, ipart+5000-1)
!
!===> ... Get standard range to be controlled
!
         irange (1, :) = max (ijk(:,i),          grid_valid_shape (1, :))
         irange (2, :) = min (ijk(:,i)+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)
!
!===> ... Set indices to be controlled
!
         if (nref == nrefd) then
            do n = 1, nref
               ijkdst (:, n) = ijk (:, i) + 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

                  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 + coords3 (i)
!
!cdir vector
         do n = 1, nref
            val = (  sin_values_lat(ijkdst (2,n)) * sin_search(i,lat)  &
                  +  cos_values_lat(ijkdst (2,n)) * cos_search(i,lat)  &
                  * (cos_values_lon(ijkdst (1,n)) * cos_search(i,lon)  &
                     +sin_values_lon(ijkdst (1,n)) * sin_search(i,lon)) )

#ifdef PRISM_ASSERTION
            if (val < acosm1 - eps .or. val > acosp1 + eps) then
               print *, i, 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 +                                 &
                           ( coords3 (i)-z_coords(ijkdst (3,n)))**2
#ifdef DEBUG2
            if (i == 2) then
               print *, 'dist, diff', dist, &
                     coords3 (i)-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 DEBUG2
         if (i == 2) then
            print *, 'psmile_neigh_near_3d_irr3_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
!
         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
!
            cycle
         endif
!
!===> ... Sort distances (silly sort)
!         Da muss noch etwas besseres her !
!
10       continue
!
         nfound = min (num_neigh, nref)
         if (nfound < num_neigh) then
            distance(nfound+1:num_neigh) = dble_huge

            do n = nfound+1, num_neigh
               ijkdst (:, n) = grid_valid_shape (1, :) - 1
            end do
         endif
!
         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
!
!===> ... Use multigrid info for search in order regions
!
#ifdef TODO
         dist_control = distance (num_neigh)
!
         call psmile_mg_search_nneigh_3d_irr3_dble (distance, ijkdst, num_neigh, &
                  grid_id, ierror)
#endif /* TODO */
!
!===> ... End of search; fill neighbors_3d
!
         do n = 1, num_neigh
            neighbors_3d (:, i, n) = ijkdst (:, n)
         end do
!
      end do ! i
      end do
!
!===> Deallocate memory
!
      Deallocate (ijk)
      Deallocate (cos_search)
      Deallocate (sin_search)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      end subroutine PSMILe_Neigh_near_3d_irr3_dble