psmile_neigh_near_3d_reg_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_reg_dble
!
! !INTERFACE:

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

      Integer, Intent (In)            :: nprev
!
!     Number of locations for which neighbours were already searched.
!
      Type (integer_vector), Intent (In) :: srclocs (ndim_3d)

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

!     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_reg_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                :: sin_search_lon (nlocs(1))
      Double Precision                :: cos_search_lon (nlocs(1))
      Double Precision                :: sin_search_lat (nlocs(2))
      Double Precision                :: cos_search_lat (nlocs(2))
!
      Integer                         :: nrefv (1:ndim_3d)
!
!     ... (1d) Indices in (entire) data to be extra searched
!
      Integer                         :: leni, lenij
      Integer                         :: np
!
!     ... 3d Indices in local data to be extra searched
!
      Integer                         :: i, j, 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
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_near_3d_reg_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_reg_dble" is designed for
! (*) grids of type "PRISM_Reglonlatvrt" and
! (*) source locations "srclocs" which are separately stored as
!     3 1-dimensional vectors
!     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_reg_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_near_3d_reg_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) * nlocs (3)) 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 ?
!
         do j = 1, ndim_3d
!cdir vector
            do i = 1, nlocs (j)
!
            if (srclocs(j)%vector(i) < grid_valid_shape(1,j) - 1 .or. &
                srclocs(j)%vector(i) > grid_valid_shape(2,j)) exit
            end do
!
         if (i <= nlocs(j)) then
            print *, "Incorrect index in j, i", j, i, srclocs(j)%vector(i)
            call psmile_assert (__FILE__, __LINE__, &
                                'wrong index')
         endif
         end do
#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.
!-----------------------------------------------------------------------
!
!===> Note: It's used that the 3d-array coords? were generated by
!           3 single regular vectors for PRISM Driver routines
!
!===> Compute values required to compute distances
!
!cdir vector
         do i = 1, nlocs (1)
         sin_search_lon (i) = coords1 (i, 1, 1) * dble_deg2rad
         end do
!
!cdir vector
         do j = 1, nlocs (2)
         sin_search_lat (j) = coords2 (1, j, 1) * dble_deg2rad
         end do
!
      cos_search_lon = cos (sin_search_lon)
      cos_search_lat = cos (sin_search_lat)
!
      sin_search_lon = sin (sin_search_lon)
      sin_search_lat = sin (sin_search_lat)
!
!-----------------------------------------------------------------------
!     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
      np = nprev
!
         do k = 1, nlocs (3)
!
!
!===>    ... Get coarse grid index on level nlev-2 and
!            standard range to be controlled
!
         if (search_mode < 3) then
            ijk(3) = srclocs(3)%vector(k)
            irange (:, 3) = ijk (3)
         else
            ijk(3) = ijk0(3) + ((srclocs(3)%vector(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))
         endif
!
         z = coords3 (1, 1, k)
         fac = dble_earth + z
!
            do j = 1, nlocs (2)
!
!===>       ... Get coarse grid index on level nlev-2 and
!               standard range to be controlled
!
            if (search_mode < 2) then
               ijk(2) = srclocs(2)%vector(j)
               irange (:, 2) = ijk (2)
            else
               ijk(2) = ijk0(2) + ((srclocs(2)%vector(j) - 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
!
               do i = 1, nlocs (1)
!
!===>          ... Get coarse grid index on level nlev-2 and
!                  standard range to be controlled
!
               ijk(1) = ijk0(1) + ((srclocs(1)%vector(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) * &
                      (irange (2,3) - irange (1,3) + 1)
!
               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 = 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
!
!cdir vector
                  do n = 1, nref
                  val = (  sin_values_lat(ijkdst (2,n)) * sin_search_lat(j)  &
                        +  cos_values_lat(ijkdst (2,n)) * cos_search_lat(j)  &
                        * (cos_values_lon(ijkdst (1,n)) * cos_search_lon(i)  &
                          +sin_values_lon(ijkdst (1,n)) * sin_search_lon(i)) )
   
#ifdef PRISM_ASSERTION
                  if (val < acosm1 - eps .or. val > acosp1 + eps) then
                     print *, i, j, 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. j == 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. j == 1 .and. k == 1) then
                  print *, 'psmile_neigh_near_3d_reg_dble: ictl =', i, j, 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_reg_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 ! j
         end do ! k
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      end subroutine PSMILe_Neigh_near_3d_reg_dble