psmile_neigh_near_irreg2_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_irreg2_dble
!
! !INTERFACE:

      subroutine psmile_neigh_near_irreg2_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, 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_irreg2_dble

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

!     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), 
                                                  coords_shape(1,2): 
                                                  coords_shape(2,2))
      Double Precision, Intent (In)   :: y_coords(coords_shape(1,1): 
                                                  coords_shape(2,1), 
                                                  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 (grid_valid_shape(1,1): 
                                                     grid_valid_shape(2,1), 
                                                     grid_valid_shape(1,2): 
                                                     grid_valid_shape(2,2), 
                                                     2)
!
!     Sin values of Longitudes and Latitudes (x_coords, y_coords)
!
      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), 
                                                     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)

      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_irreg2_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
!
!  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              :: 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, iz, n, np
      Double Precision                :: z
!
!     ... For locations controlled
!
      Double Precision, Pointer       :: sin_search (:, :)
      Double Precision, Pointer       :: cos_search (:, :)
!
      Integer                         :: ijkctl (ndim_3d, nref_3d)
      Integer                         :: nrefv (2:ndim_3d)
      Integer                         :: ijkdst (ndim_3d, nref_3d)
      Integer                         :: irange (2, ndim_3d)
      Integer                         :: itemp (ndim_3d)
!
      Logical                         :: mask_dist (nref_3d)
      Double Precision                :: distance (nref_3d), dist, fac, val
      Double Precision                :: dist_control
!
      Integer                         :: ijk (ndim_3d), ijk0 (ndim_3d)
      Integer                         :: ii, jj, kk, leni, lenij, lenk
      Integer                         :: nref, nrefd
!
      Integer                         :: nfound
      Integer                         :: nmin (1)
!
!     ...
!
      Double Precision                :: dble_huge
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_near_irreg2_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 2 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).
!
! Subroutine "PSMILe_Neigh_near_irreg2_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_irreg2_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_near_irreg2_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
!  Relative control indices 
!      
! Schoener waere eine liste, die weniger sortierung nachher verlangt
! (dh. die nach den Summe der relativen Indices sortiert ist).
! Vielleicht 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_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 < 2 .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) - 1 .or. &
             srcloc(1,i) > grid_valid_shape(2,1)     .or. &
             srcloc(2,i) < grid_valid_shape(1,2) - 1 .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 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(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) * dble_deg2rad
         end do
!
!cdir vector
         do i = 1, nlocs (1)
         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, &
!
!-----------------------------------------------------------------------
!
      lenk = 1
      ijk0 = Grids(grid_id)%ijk0
      np = nprev
!
         do iz = 1, nlocs (2)
!
!===> ... Get coarse grid index on level nlev-2
!
         if (search_mode == 2) then
            ijk(3) = srclocz(iz)
            irange (:, 3) = srclocz(iz)
         else
!
            ijk(3) = ijk0(3) + ((srclocz(iz) - 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 iz
!       berechnen, falls mask = true oder nicht available
!       Vorsicht bei Reihenfolge np+i
!
         z = coords3 ((iz-1)*nlocs (1) + 1)
         fac = dble_earth + z
!
            do i = 1, nlocs (1)
!
!===>       ... Get coarse grid index on level nlev-2 and
!               standard range to be controlled
!               Get standard range to be controlled
!               Note: irange (:, 3) is already set above for search_mode == 2
!
            ijk(1:2) = ijk0(1:2) + ((srcloc(:,i) - ijk0(1:2)) / 4) * 4
!
            irange (1, 1:2) = max (ijk(1:2),   grid_valid_shape (1, 1:2))
            irange (2, 1:2) = min (ijk(1:2)+3, grid_valid_shape (2, 1:2))
!
            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(ijkdst (1,n),                              &
                                     ijkdst (2,n), lat) * sin_search(i, lat)    &
                       +  cos_values(ijkdst (1,n),                              &
                                     ijkdst (2,n), lat) * cos_search(i, lat)    &
                       * (cos_values(ijkdst (1,n),                              &
                                     ijkdst (2,n), lon) * cos_search(i, lon)    &
                         +sin_values(ijkdst (1,n),                              &
                                     ijkdst (2,n), lon) * sin_search(i, lon)) )

#ifdef PRISM_ASSERTION
               if (val < acosm1 - eps .or. val > acosp1 + eps) then
                  print *, i, iz, 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. iz == 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. iz == 1) then
               print *, 'psmile_neigh_near_irreg2_dble: ictl =', i, iz, 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
!
            dist_control = distance (num_neigh)
!
#ifdef TODO
            call psmile_mg_search_nneigh_irreg2_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 ! iz
!
!===> 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_irreg2_dble')
9980 format (1x, a, ': psmile_neigh_near_irreg2_dble: eof, ierror =', i3)

      end subroutine PSMILe_Neigh_near_irreg2_dble