psmile_neigh_near_3d_irr2_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_near_3d_irr2_dble
!
! !INTERFACE:

      subroutine psmile_neigh_near_3d_irr2_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, srclocz, nlocs,         &
                           nloc, nprev,                    &
                           neighbors_3d, num_neigh,        &
                           extra_search, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_near_3d_irr2_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)            :: nlocs (ndim_2d)
!
!     Number of locations to be transferred with
!     nlocs (1) = Number of locations in lon lat direction
!     nlocs (2) = Number of locations in z       direction

      Integer, Intent (In)            :: nprev
!
!     Number of locations for which neighbours were already searched.
!
      Integer, Intent (In)            :: srcloc  (ndim_2d, nlocs (1))
      Integer, Intent (In)            :: srclocz (nlocs (2))

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

!     Coordinates of the target points (which were searched and found)
!     Note: Target grids are full arrays which were generated 
!           out of partial arrays.

      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_irr2_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                         :: leni, lenij, lenk
      Integer                         :: np
!
!     ... 3d Indices in local data to be extra searched
!
      Integer                         :: i, k
      Integer                         :: ii, jj, kk
      Integer                         :: n, nref, nrefd
!
      Integer                         :: ijk (ndim_3d)
      Integer                         :: ijk0 (ndim_3d)
      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)
!
!     ... for nearest neighbour control
!
      Integer                         :: nfound
      Integer                         :: nmin (1)
!
      Double Precision                :: dble_huge
      Double Precision                :: dist, fac, val, z
      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_irr2_dble" searches the next "num_neigh"
! nearest neighbours on the method-grid (x_coords, y_coords, z_coords)
! on the 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_irr2_dble" is designed for
! (*) grids of type "PRISM_Reglonlatvrt" and
! (*) source locations "srcloc, srclocz" which are stored as
!     a 2-dimensional vector for the lon lat directions and 
!     a 1-dimensional vector for the z-direction
!     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_irr2_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_near_3d_irr2_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
!  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
!
#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 < nprev + nlocs (1) * nlocs (2)) then
         print *, trim(ch_id), " nloc, nprev, nlocs ", nloc, nprev, nlocs
         call psmile_assert (__FILE__, __LINE__, &
                             'nloc < nprev + PRODUCT (nlocs) ')
      endif
!
!===> Are the locations within the correct shape ?
!
!cdir vector
         do i = 1, nlocs (1)
!
         if (srcloc(1,i) < grid_valid_shape(1,1) .or. &
             srcloc(1,i) > grid_valid_shape(2,1) .or. &
             srcloc(2,i) < grid_valid_shape(1,2) .or. &
             srcloc(2,i) > grid_valid_shape(2,2)) exit
         end do
!
      if (i <= nlocs(1)) then
         print *, "Incorrect index in srcloc, i", i, srcloc(:, i)
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index in srcloc')
      endif

!cdir vector
         do i = 1, nlocs (2)
         if (srclocz(i) < grid_valid_shape(1,3) .or. &
             srclocz(i) > grid_valid_shape(2,3)) exit
         end do
!
      if (i <= nlocs(2)) then
         print *, "Incorrect index in srclocz(2), i", i, srclocz(i)
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index in srclocz')
      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 points 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(nlocs(1), lat), STAT = ierror)

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

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

      Allocate (cos_search(nlocs(1), lat), STAT = ierror)

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

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'cos_search', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!===> Compute values required to compute distances
!     Note: Target grids are full arrays which were generated 
!           out of partial arrays.
!
!cdir vector
         do i = 1, nlocs (1)
         sin_search (i, lon) = coords1 (i,1) * dble_deg2rad
         end do
!
!cdir vector
         do i = 1, nlocs (1)
         sin_search (i, lat) = coords2 (i,1) * 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, &
!
!-----------------------------------------------------------------------
!
      lenk = 1
      ijk0 = Grids(grid_id)%ijk0
      np = nprev
!
         do k = 1, nlocs (2)
!
!===>    ... Get coarse grid index on level nlev-2 and
!            standard range to be controlled
!
         if (search_mode < 3) then
            ijk(3) = srclocz(k)
            irange (:, 3) = ijk (3)
         else
            ijk(3) = ijk0(3) + ((srclocz(k) - ijk0(3)) / 4) * 4
            irange (1, 3) = max (ijk (3),   grid_valid_shape (1, 3))
            irange (2, 3) = min (ijk (3)+3, grid_valid_shape (2, 3))
            lenk = irange (2,3) - irange (1,3) + 1
         endif
!
! TODO: Loops umdrehen und die sin und cos values nur einmal fuer k
!       berechnen, falls mask = true oder nicht available
!       Vorsicht bei Reihenfolge np+i
!
         z = coords3 (1, k)
         fac = dble_earth + z
!
            do i = 1, nlocs (1)
!
!===>       ... Get coarse grid index on level nlev-2 and
!               standard range to be controlled
!
            if (search_mode < 2) then
               ijk(2) = srcloc(2,i)
               irange (:, 2) = ijk (2)
            else
               ijk(2) = ijk0(2) + ((srcloc(2,i) - ijk0(2)) / 4) * 4
               irange (1, 2) = max (ijk (2),   grid_valid_shape (1, 2))
               irange (2, 2) = min (ijk (2)+3, grid_valid_shape (2, 2))
            endif
!
            ijk(1) = ijk0(1) + ((srcloc(1,i) - ijk0(1)) / 4) * 4
!
            irange (1, 1) = max (ijk(1),   grid_valid_shape (1, 1))
            irange (2, 1) = min (ijk(1)+3, grid_valid_shape (2, 1))
!
            nref = (irange (2,1) - irange (1,1) + 1) * &
                   (irange (2,2) - irange (1,2) + 1) * &
                   lenk
!
!===> ... Set indices to be controlled
!
            if (nref == nrefd) then
                  do n = 1, nrefd
                  ijkdst (:, n) = ijk + ijkctl (:, n)
                  end do
            else
! 
               leni  =  irange (2, 1) - irange (1, 1) + 1
               lenij = (irange (2, 2) - irange (1, 2) + 1) * leni 
!
               do kk = 0, lenk-1
               n = kk*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*lenij+1:(kk+1)*lenij) = irange (1,3) + kk
               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
!
!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, k, 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-z_coords(ijkdst (3,n)))**2
#ifdef DEBUG2
               if (i == 1 .and. k == 1) then
                  print *, 'dist, diff', dist, &
                           z-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 == 1 .and. k == 1) then
               print *, 'psmile_neigh_near_irreg2_dble: ictl =', i, k, np
               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 (:, np+i, 1) = ijkdst (:, nmin(1))
!
                  do n = 2, num_neigh
                  neighbors_3d (:, np+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
!
#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_irr2_dble (distance, ijkdst, num_neigh, &
                        grid_id, ierror)
#endif /* TODO */
!
!===> ... End of search; fill neighbors_3d
!
               do n = 1, num_neigh
               neighbors_3d (:, np+i, n) = ijkdst (:, n)
               end do
!
            end do ! i

         np = np + nlocs (1)
         end do ! k
!
!===> Deallocate memory
!
      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_irr2_dble')
9980 format (1x, a, ': psmile_neigh_near_3d_irr2_dble: eof, ierror =', i3)

      end subroutine PSMILe_Neigh_near_3d_irr2_dble