psmile_neigh_nearx_3d_reg_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_nearx_3d_reg_dble
!
! !INTERFACE:

      subroutine psmile_neigh_nearx_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_nearx_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)
!
!     Total  Number of locations to be transferred

      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_nearx_3d_reg_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
!  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
!
      Integer, Parameter              :: lon = 1
      Integer, Parameter              :: lat = 2
!
! !LOCAL VARIABLES
!
!     ... For locations controlled
!
      Double Precision, Pointer       :: sin_search (:, :)
      Double Precision, Pointer       :: cos_search (:, :)
      Double Precision, Pointer       :: z_search (:)
!
!     ... (1d) Indices in (entire) data to be extra searched
!
#ifdef PRISM_ASSERTION
      Integer                         :: i
#endif
      Integer                         :: j
      Integer                         :: ibeg, iend
      Integer                         :: jbeg, jend
      Integer                         :: lenij
!
      Integer, Pointer                :: indices (:)
      Integer                         :: nprev1
!
!     ... 3d Indices in local data to be extra searched
!
      Integer                         :: ijk0 (ndim_3d)
      Integer, Pointer                :: ijk (:, :)
      Integer                         :: width (ndim_3d)
!
!     ... for error handling
!
      Integer, Parameter              :: nerrp = 2
      Integer                         :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Neigh_nearx_3d_reg_dble" searches the next "num_neigh"
! nearest neighbours on the method-grid (x_coords, y_coords, z_coords)
! for the extra indices specified in "extra_search%indices".
!
! 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_nearx_3d_reg_dble" is designed
! (*) for source grids of type "PRISM_Reglonlatvrt" and
! (*) source locations "srcloc" which are stored as a 3 single 1d vectors
!     containing the 3d indices of the cells found for the target points
!     "coords" to be searched.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_neigh_nearx_3d_reg_dble.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_neigh_nearx_3d_reg_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
!                            
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      ierror = 0
!
      nprev1 = nprev + 1
!
      ibeg = nprev + 1
      iend = nprev + nlocs(1)*nlocs(2)*nlocs(3)
!
      indices => extra_search%indices
!
!===> Get range in indices to be controlled
!     Note: indices are locations in the full 3d cube and
!           contains the previous numbers
!           TODO: DAS SOLLTE in lokale Indices geaendert werden
!
         do jbeg = 1, extra_search%n_extra
         if (indices(jbeg) >= ibeg) exit
         end do
!
      if (jbeg > extra_search%n_extra) return
!
         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
!
      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 ?
!     Note: the scrloc's are source locations in the
!           method grid which has an virtual cell 0.
!
         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 in srcloc')
         endif
         end do
#endif
!
!-----------------------------------------------------------------------
!     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(jbeg:jend, lat), &
                cos_search(jbeg:jend, lat), &
                  z_search(jbeg:jend),      STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = (jend-jbeg+1) * (lat * 2 + 1)

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'sin_search/cos_search/z_search', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
      Allocate (ijk(ndim_3d, jbeg:jend), STAT = ierror)

      if ( ierror > 0 ) then
         ierrp (1) = ierror
         ierrp (2) = (jend-jbeg+1) * ndim_3d

         ierror = PRISM_Error_Alloc
         call psmile_error ( ierror, 'ijk', &
                             ierrp, 2, __FILE__, __LINE__ )
         return
      endif
!
!===> Compute 3d Indices
!
      lenij = nlocs(1) * nlocs (2)
!
      ijk (3, jbeg:jend) =    (indices(jbeg:jend)-nprev1)/lenij + 1
      ijk (2, jbeg:jend) =   ((indices(jbeg:jend)-nprev1) - &
                                (ijk(3,jbeg:jend)-1)*lenij)/nlocs(1) + 1
      ijk (1, jbeg:jend) = mod(indices(jbeg:jend)-nprev1, nlocs(1)) + 1

#ifdef PRISM_ASSERTION
!cdir vector
         do j = jbeg, jend
         if (ijk(1,j) < 1 .or. ijk(1,j) > nlocs(1) .or. &
             ijk(2,j) < 1 .or. ijk(2,j) > nlocs(2) .or. &
             ijk(3,j) < 1 .or. ijk(3,j) > nlocs(3)) exit
         end do

      if (j <= jend) then
         print *, "Incorrect index in ijk", j, jbeg, jend, ijk(:, j)
         call psmile_assert (__FILE__, __LINE__, &
                                'wrong index in ijk')
      endif
#endif
!
!===> Compute values required to compute distances
!     Note: Target grids are full arrays which were generated 
!           out of partial arrays.
!
!cdir vector
         do j = jbeg, jend
         sin_search (j, lon) = coords1 (ijk(1,j), ijk(2,j), ijk(3,j)) &
                             * dble_deg2rad
!        sin_search (j, lon) = coords1 (ijk(1,j), 1, 1) * dble_deg2rad
         end do
!
!cdir vector
         do j = jbeg, jend
         sin_search (j, lat) = coords2 (ijk(1,j), ijk(2,j), ijk(3,j)) &
                             * dble_deg2rad
!        sin_search (j, lon) = coords1 (1, ijk(2,j), 1) * dble_deg2rad
         end do
!
      cos_search = cos (sin_search)
      sin_search = sin (sin_search)
!
!cdir vector
         do j = jbeg, jend
         z_search (j) = coords3 (ijk(1,j), ijk(2,j), ijk(3,j))
!        z_search (j) = coords3 (1, 1, ijk(3,j))
         end do
!
!-----------------------------------------------------------------------
!     Compute coarse grid indices
!-----------------------------------------------------------------------
!
      ijk0 = Grids(grid_id)%ijk0
!
!===> Get coarse grid index on level nlev-2.
!     By coarsening with a faktor of 4
!
      width (1:search_mode) = 3
      width (search_mode+1:ndim_3d) = 0
!
      ijk (1, jbeg:jend) = ijk0(1) + ((srclocs(1)%vector(ijk(1, jbeg:jend)) - &
                                       ijk0(1)) / 4) * 4
!
      if (search_mode >= 2) then
         ijk (2, jbeg:jend) = ijk0(2) + ((srclocs(2)%vector(ijk(2, jbeg:jend)) - &
                                         ijk0(2)) / 4) * 4
      else
         ijk (2, jbeg:jend) = srclocs(2)%vector(ijk(2, jbeg:jend))
      endif
!
      if (search_mode >= 3) then
         ijk (3, jbeg:jend) = ijk0(3) + ((srclocs(3)%vector(ijk(3, jbeg:jend)) - &
                                         ijk0(3)) / 4) * 4
      else
         ijk (3, jbeg:jend) = srclocs(3)%vector(ijk(3, jbeg:jend))
      endif
!
!===> Compute nearest neigbours for all points
!     specified in "indices(jbeg:jend)"
!
      call psmile_neigh_nearx_sub_reg_dble (               &
                           grid_id,                        &
                           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,  &
                           neighbors_3d, num_neigh, nloc,  &
                           extra_search,                   &
                           ijk, sin_search, cos_search,    &
                           z_search, jbeg, jend, width, ierror)
!
!===> Deallocate memory
!
      Deallocate (ijk)
      Deallocate (cos_search, sin_search, z_search)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      end subroutine PSMILe_Neigh_nearx_3d_reg_dble