psmile_neigh_nearestx_3d_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_nearestx_3d_dble
!
! !INTERFACE:

      subroutine psmile_neigh_nearestx_3d_dble (           &
                           grid_id,                        &
                           coords1, coords2, coords3,      &
                           x_coords, y_coords, z_coords,   &
                           coords_shape,                   &
                           mask_array, mask_shape,  mask_available,    &
                           sin_values, cos_values,         &
                           grid_valid_shape,               &
                           srcloc, nloc, nprev, nsearch,   &
                           neighbors_3d, num_neigh,        &
                           extra_search, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Neigh_nearestx_3d_dble

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

!     Link to the info on the component in which the donor
!     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
!     run.
!
      Integer, Intent (In)            :: srcloc (ndim_3d, nloc)
!
!     Locations to be transferred

!     Indices of the grid cell in which the point was found.
!     The indices are relative to shape.

      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), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))
      Double Precision, Intent (In)   :: y_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))
      Double Precision, Intent (In)   :: z_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2), 
                                                  coords_shape(1,3): 
                                                  coords_shape(2,3))

!     Coordinates of the method

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

!     Dimension of (method) mask array

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

      Logical, Intent (In)            :: mask_available

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

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

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

      Double Precision, Intent (In)   :: sin_values (grid_valid_shape(1,1): 
                                                     grid_valid_shape(2,1), 
                                                     grid_valid_shape(1,2): 
                                                     grid_valid_shape(2,2), 
                                                     grid_valid_shape(1,3): 
                                                     grid_valid_shape(2,3), 
                                                     2)
!
!     Sin values of Longitudes and Latitudes (x_coords, y_coords)
!
      Double Precision, Intent (In)   :: cos_values (grid_valid_shape(1,1): 
                                                     grid_valid_shape(2,1), 
                                                     grid_valid_shape(1,2): 
                                                     grid_valid_shape(2,2), 
                                                     grid_valid_shape(1,3): 
                                                     grid_valid_shape(2,3), 
                                                     2)
!
!     Cos values of Longitudes and Latitudes (x_coords, y_coords)
!

      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_nearestx_3d_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NREFD = Number of locations to be controlled
!        = 4 ** ndim_3d
!
!  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              :: nrefd = 4 * 4 * 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 searched
!
      Integer                         :: i, ibeg, iend, n
      Integer                         :: j, jbeg, jend
!
      Integer, Pointer                :: indices (:)
!
!     ... For locations controlled
!
      Double Precision, Pointer       :: sin_search (:, :)
      Double Precision, Pointer       :: cos_search (:, :)
      Double Precision, Pointer       :: z_search (:)
!
      Integer                         :: ijkctl (ndim_3d, nrefd)
      Integer                         :: ijkdst (ndim_3d, nrefd)
      Integer                         :: irange (2, ndim_3d)
      Integer                         :: itemp (ndim_3d)
!
      Logical                         :: mask_dist (nrefd)
      Double Precision                :: distance (nrefd), dist, fac, val
      Double Precision                :: dist_control
!
      Integer                         :: ijk (ndim_3d), ijk0 (ndim_3d)
      Integer                         :: leni, lenij
      Integer                         :: ii, jj, kk, nref
!
      Integer                         :: nfound
      Integer                         :: nmin (1)
!
!     ... 
!
      Double Precision                :: dble_huge
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_nearestx_3d_dble" searches the next "num_neigh"
! neighbours on the method-grid (x_coords, y_coords, z_coords)
! for the extra indices specified in "extra_search%indices".
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_nearestx_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_nearestx_3d_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,nrefd) &
                            /  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,nrefd) &
            /  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/

!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      ibeg = nprev + 1
      iend = nprev + nsearch
!
      indices => extra_search%indices
!
!===> Get range in "indices(:)" to be controlled
!
!cdir vector
         do jbeg = 1, extra_search%n_extra
         if (indices(jbeg) >= ibeg) exit
         end do
!
      if (jbeg > extra_search%n_extra) return
!
!cdir vector
         do jend = jbeg, extra_search%n_extra
         if (indices(jend) > iend) exit
         end do
      jend = jend - 1
!
      if (jbeg > jend) return
!
      ibeg = indices (jbeg)
      iend = indices (jend)
!
#ifdef PRISM_ASSERTION
!
!===> Internal control
!
!cdir vector
         do j = jbeg, jend
         i = indices(j)
!
         if (srcloc(1,i) < coords_shape(1,1) .or. &
             srcloc(1,i) > coords_shape(2,1) .or. &
             srcloc(2,i) < coords_shape(1,2) .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 (j <= jend) then
         print *, i, srcloc(:, i), coords_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index')
      endif
#endif
!
      dble_huge = huge (distance(1))
!
!-----------------------------------------------------------------------
!     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(jbeg:jend, 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(jbeg:jend, 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 (z_search(jbeg:jend), STAT = ierror)

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

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'z_search', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!cdir vector
         do j = jbeg, jend
         sin_search (j, lon) = coords1 (indices(j)) * dble_deg2rad
         end do
!
!cdir vector
         do j = jbeg, jend
         sin_search (j, lat) = coords2 (indices(j)) * dble_deg2rad
         end do
!
      cos_search = cos (sin_search)
      sin_search = sin (sin_search)
!
!cdir vector
         do j = jbeg, jend
         z_search (j) = coords3 (indices(j))
         end do
!
!-----------------------------------------------------------------------
!     Compute distances
!
!     dist = (dble_earth + z_search)*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 + (z_search-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
!
         do j = jbeg, jend
         i = indices(j)
!
!===> Get coarse grid index on level nlev-2
!
         ijk = ijk0 + ((srcloc(:,i) - ijk0) / 4) * 4
!
!===> Get standard range to be controlled
!
         irange (1, :) = max (ijk,   grid_valid_shape (1, :))
         irange (2, :) = min (ijk+3, 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)
!
         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
!
         fac = dble_earth + z_search (j)
!cdir vector
            do n = 1, nref
            val   = (  sin_values(ijkdst (1,n),                               &
                                  ijkdst (2,n),                               &
                                  ijkdst (3,n), lat) * sin_search(j, lat)     &
                    +  cos_values(ijkdst (1,n),                               &
                                  ijkdst (2,n),                               &
                                  ijkdst (3,n), lat) * cos_search(j, lat)     &
                    * (cos_values(ijkdst (1,n),                               &
                                  ijkdst (2,n),                               &
                                  ijkdst (3,n), lon) * cos_search(j, lon)     &
                      +sin_values(ijkdst (1,n),                               &
                                  ijkdst (2,n),                               &
                                  ijkdst (3,n), lon) * sin_search(j, lon)) )
#ifdef PRISM_ASSERTION
            if (val < acosm1 - eps .or. val > acosp1 + eps) then
               print *, i, j, val
               call psmile_assert (__FILE__, __LINE__, &
                                   'invalid value for acos')
            endif
#endif
            val = min (acosp1, val)
            val = max (acosm1, val)

            dist = fac * acos (val)

            distance (n) = dist*dist +                                 &
                           ( z_search (j)-z_coords(ijkdst (1,n),       &
                                                   ijkdst (2,n),       &
                                                   ijkdst (3,n)))**2
#ifdef DEBUG2
         if (i == 2) then
            print *, 'dist, diff', dist, &
               z_search (j)-z_coords(ijkdst (1,n), ijkdst (2,n), 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_nearestx_3d_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
!
#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
!
         dist_control = distance (num_neigh)
!
#ifdef TODO
         call psmile_mg_search_nneigh_3d_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 ! j
!
!===> Deallocate memory
!
      Deallocate (z_search)
      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_nearestx_3d_dble')
9980 format (1x, a, ': psmile_neigh_nearestx_3d_dble: eof, ierror =', i3)

      end subroutine PSMILe_Neigh_nearestx_3d_dble