psmile_mg_search.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2006-2010, NEC Europe Ltd., London, UK.
! Copyright 2011, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: psmile_mg_search
!
! !INTERFACE:

      subroutine psmile_mg_search (comp_info, grid_id, search_data, &
                                   requires_conserv_remap, found,   &
                                   locations, len, tol, ierror)
!
! !USES:
!
      use PRISM_constants
      use PSMILe, dummy_interface => psmile_mg_search
      use psmile_grid, only : get_num_independent_dims

      implicit none
!
! !INPUT PARAMETERS:
!
      Type (Enddef_comp), Intent (In)           :: comp_info
!
!     Information on the component
!
      Integer, Intent (In)                      :: grid_id
!
!     Local grid id
!
      Double Precision, Intent (In)             :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0
!
      Logical, Intent (In)                      :: requires_conserv_remap
!
!     Are the results later used for conservative remapping (only
!     significant for reduced gaussian grids)
!
!
! !INPUT/OUTPUT PARAMETERS:
!
      Type (Enddef_search_data), Intent (InOut) :: search_data
!
!     Data on the points to be searched
!
! !OUTPUT PARAMETERS:
!
      Type (integer_vector), Intent (Out)       :: found     (search_data%npart, ndim_3d)
      Type (integer_vector), Intent (Out)       :: locations (search_data%npart, ndim_3d)
                                               
      Integer, Intent (Out)                     :: len  (search_data%npart, ndim_3d)
                                               
      Integer, Intent (Out)                     :: ierror

!     Returns the error code of psmile_mg_search;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
      Integer                     :: i
      Real                        :: rtol
#if defined ( PRISM_QUAD_TYPE )
      Real (kind=PRISM_QUAD_TYPE) :: qtol
#endif
!
!     ... for grids and components
!
      Integer                     :: nlev1
      Integer                     :: datatype
!
!     ... for intersections
!
      Integer                     :: ipart, npart
!
!     ... locations searched
!
      Integer                     :: n_vec
!
!     ... for error parameters
!
      Integer, parameter          :: nerrp = 2
      Integer                     :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "psmile_mg_search" creates the found and locations arrays and filles them
! with the inital results of the multigrid search
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.21    H. Ritzdorf  created as psmile_search_donor_cells
! 03.05.2011  M. Hanke     created by moving parts of psmile_search_donor_cells
!                          to this new files/routine
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_mg_search.F90 3215 2011-06-07 10:00:41Z coquart $
! $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_search.F90 3215 2011-06-07 10:00:41Z coquart $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)
      call psmile_flushstd
#endif /* VERBOSE */

      rtol = tol

#if defined ( PRISM_QUAD_TYPE )
      qtol = tol
#endif

      datatype  = Grids(grid_id)%corner_pointer%corner_datatype

      npart = search_data%npart
      nlev1 = - (Grids(grid_id)%nlev + 1) ! Outside

      ! for reduced gauss grids an auxiliary grid of type PRISM_Reglonlatvrt
      ! is used for the search
      if (Grids(grid_id)%grid_type == PRISM_Gaussreduced_regvrt) then
         n_vec = get_num_independent_dims(PRISM_Reglonlatvrt)
      else
         n_vec = get_num_independent_dims(Grids(grid_id)%grid_type)
      endif

! -----------------------------------------------------------------------
!      calculate the size of the found and location arrays
!      and allocate them depending on number of independent
!      of the local grid type
! -----------------------------------------------------------------------

      select case ( n_vec )

      case (ndim_3d) ! 3 independent dimension
                     ! regular in all 3 directions

!    len (ipart, i) = Length of last dimension of vectors
!                     found (ipart,i) and locations (ipart,i)
!                   = Number of coordinates to be searched in each direction
!


         ! check whether the remote grid is among the supported grid types
         if ( .not. any (search_data%grid_type == (/PRISM_Reglonlatvrt,     &
                                                    PRISM_Irrlonlat_regvrt, &
                                                    PRISM_Irrlonlatvrt,     &
                                                    PRISM_Gaussreduced_regvrt/))) then
            ierrp (1) = search_data%grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'unsupported grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
         endif

         ! compute the size of the found and location arrays
         ! 1 vector for each direction of the source grid
         do ipart = 1, npart
            len (ipart,:) = search_data%dim_size(:,ipart)
         enddo

         ! For reduced gauss grids in combination with conservative remapping we have
         ! two points per cell (lower left and upper right).
         ! In order for the multigrid search to take all points into account we need
         ! to extent the search range.
         if ( search_data%grid_type == PRISM_Gaussreduced_regvrt .and. &
              requires_conserv_remap) then

! TODO: Isn't the doubling already done in dim_size?
            len(:,1) = 2 * len(:,1)
            len(:,2) = len(:,1)

            search_data%range(2, 1, 1:npart) = 2 * search_data%range(2, 1, 1:npart) - &
                                                   search_data%range(1, 1, 1:npart) + 1
            search_data%shape(2, 1, 1:npart) = search_data%range(2, 1, 1:npart)
         endif

         ! allocate the found and location arrays
         do ipart = 1, npart

            do i = 1, ndim_3d
               Allocate (found(ipart,i)%vector(len(ipart,i)), &
                         locations(ipart,i)%vector(len(ipart,i)), STAT = ierror)
               if ( ierror > 0 ) then
                  ierrp (1) = ierror
                  ierrp (2) = 2*len (ipart, i)
                  ierror = PRISM_Error_Alloc
                  call psmile_error ( ierror, 'found(ipart,i)%vector, locations(ipart,i)%vector', &
                                      ierrp, 2, __FILE__, __LINE__ )
                  return
               endif
            enddo
         end do ! ipart

         ! initalise found and locations arrays
         do ipart = 1, npart
            do i = 1, ndim_3d
               found    (ipart,i)%vector (:) = nlev1
               locations(ipart,i)%vector (:) = 0
            end do
         end do ! ipart

      case (ndim_2d) ! 2 independent dimesions
                     ! Irregular in lonlat   direction
                     ! Regular   in vertical direction

!    len (ipart, i) = Length of last dimension of vectors
!                     found (ipart,i) and locations (ipart,i)
!                   = Number of coordinates to be searched with
!                     len (ipart, 1) = Number in lonlat direction
!                     len (ipart, 2) = Number in z      direction

         ! compute the size of the found and location arrays depending on the
         ! number of independent dimension of the target grid type
         ! 2 vectors for lonlat and z direction of the source grid
         if (get_num_independent_dims(search_data%grid_type) == ndim_3d) then

            do ipart = 1, npart
               len(ipart,1) = product (search_data%dim_size(1:2,ipart))
               len(ipart,2) = search_data%dim_size(3, ipart)
            enddo

         else

            do ipart = 1, npart
               len(ipart,1) = search_data%dim_size(1, ipart)
               len(ipart,2) = search_data%dim_size(3, ipart)
            end do

         endif

         ! For reduced gauss grids in combination with conservative remapping we have
         ! two points per cell (lower left and upper right).
         ! In order for the multigrid search to take all points into account we need
         ! to extent the search range.
         if ( search_data%grid_type == PRISM_Gaussreduced_regvrt .and. &
              requires_conserv_remap) then

            search_data%range(2, 1, 1:npart) = 2 * search_data%range(2, 1, 1:npart) - &
                                                   search_data%range(1, 1, 1:npart) + 1
            search_data%shape(2, 1, 1:npart) = search_data%range(2, 1, 1:npart)
         endif

         do ipart = 1, npart
            Allocate (found(ipart,1)%vector(len(ipart,1)), &
                      locations(ipart,1)%vector(ndim_2d*len(ipart,1)), &
                      found(ipart,2)%vector(len(ipart,2)), &
                      locations(ipart,2)%vector(len(ipart,2)), STAT = ierror)
            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = (ndim_2d+1) * len(ipart,1) + 2*len (ipart,2)
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'found(ipart,1:2)%vector, locations(ipart,1:2)%vector', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
            Nullify (found(ipart,3)%vector, locations(ipart,3)%vector)
         end do ! ipart

         do ipart = 1, npart
            found    (ipart,1)%vector (:) = nlev1
            found    (ipart,2)%vector (:) = nlev1

            locations(ipart,1)%vector (:) = 0
            locations(ipart,2)%vector (:) = 0
         end do ! ipart

      case (ndim_1d) ! 1 independent dimension
                     ! Irregular in lonlat and vertical direction

!    len (ipart, i) = Length of last dimension of vectors
!                     found (ipart,i) and locations (ipart,i)
!                   = Number of coordinates to be searched

          ierrp (1) = Grids(grid_id)%grid_type
          ierror = PRISM_Error_Internal

          call psmile_error ( ierror, 'untested part of the code', &
                              ierrp, 1, __FILE__, __LINE__ )


         ! compute the size of the found and location arrays
         ! 1 vector for all directions
         do ipart = 1, npart
            len(ipart,1) = product(search_data%range (2,1:ndim_3d,ipart) - &
                                   search_data%range (1,1:ndim_3d,ipart) + 1)
         end do

         do ipart = 1, npart
            Allocate (found(ipart,1)%vector(len(ipart,1)), &
                      locations(ipart,1)%vector(ndim_3d*len(ipart,1)), STAT = ierror)
            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = (ndim_3d+1)*len(ipart,1)
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'found(ipart,1)%vector, locations(ipart,1)%vector', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
            Nullify (found(ipart,2)%vector, locations(ipart,2)%vector, &
                     found(ipart,3)%vector, locations(ipart,3)%vector)
         end do ! ipart

         do ipart = 1, npart
            found    (ipart,1)%vector (:) = nlev1
            locations(ipart,1)%vector (:) = 0
         end do ! ipart

      case DEFAULT
!
          ierrp (1) = Grids(grid_id)%grid_type
          ierror = PRISM_Error_Internal

          call psmile_error ( ierror, 'unsupported grid type', &
                              ierrp, 1, __FILE__, __LINE__ )

      end select ! n_vec

! -----------------------------------------------------------------------
!      Do the search on the multigrid and afterwards the actual search
!      + everything else there is to do
! -----------------------------------------------------------------------
      
      select case ( n_vec )

      case (ndim_3d) ! 3 independent dimension
                     ! regular in all 3 directions

         ! Search locations separately in all 3 directions
         if (datatype == MPI_REAL) then

            do ipart = 1, npart
               do i = 1, ndim_3d
                  call psmile_mg_search_1d_real (grid_id, i,            &
                              found(ipart,i)%vector,                    &
                              locations(ipart,i)%vector,                &
                              search_data%search_real(i, ipart)%vector, &
                              search_data%dim_size(i, ipart), rtol, ierror)
                  if (ierror > 0) return
               end do
            end do ! ipart

         else if (datatype == MPI_DOUBLE_PRECISION) then

            do ipart = 1, npart
               do i = 1, ndim_3d
                  call psmile_mg_search_1d_dble (grid_id, i,            &
                              found(ipart,i)%vector,                    &
                              locations(ipart,i)%vector,                &
                              search_data%search_dble(i, ipart)%vector, &
                              search_data%dim_size(i, ipart), tol, ierror)
                  if (ierror > 0) return
               end do
            end do ! ipart

#if defined ( PRISM_QUAD_TYPE )
         else if (datatype == MPI_REAL16) then

            do ipart = 1, npart
               do i = 1, ndim_3d
                  call psmile_mg_search_1d_quad (grid_id, i,            &
                              found(ipart,i)%vector,                    &
                              locations(ipart,i)%vector,                &
                              search_data%search_quad(i, ipart)%vector, &
                              search_data%dim_size(i, ipart), qtol, ierror)
                  if (ierror > 0) return
               end do
            end do ! ipart

#endif
         endif

      case (ndim_2d) !    2 independent dimesions
                     !    Irregular in lonlat   direction
                     !    Regular   in vertical direction

         if (datatype == MPI_REAL) then

            do ipart = 1, npart
               call psmile_mg_search_2d_real (grid_id,                       &
                           found(ipart,1)%vector, locations(ipart,1)%vector, &
                           len  (ipart,1), search_data, ipart, rtol, ierror)
               if (ierror > 0) return

               call psmile_mg_search_1d_real (grid_id, 3,            &
                           found    (ipart,2)%vector,                &
                           locations(ipart,2)%vector,                &
                           search_data%search_real(3, ipart)%vector, &
                           search_data%dim_size(3, ipart), rtol, ierror)
               if (ierror > 0) return
            end do

         else if (datatype == MPI_DOUBLE_PRECISION) then

            do ipart = 1, npart

#ifdef DEBUG
           PRINT*, TRIM(ch_id), ' intersection :',ipart
           call psmile_flushstd
#endif
               call psmile_mg_search_2d_dble (grid_id,                       &
                           found(ipart,1)%vector, locations(ipart,1)%vector, &
                           len  (ipart,1), search_data, ipart, tol, ierror)
               if (ierror > 0) return

               call psmile_mg_search_1d_dble (grid_id, 3,                    &
                           found(ipart,2)%vector, locations(ipart,2)%vector, &
                           search_data%search_dble(3, ipart)%vector,         &
                           search_data%dim_size(3, ipart), tol, ierror)
               if (ierror > 0) return
            end do

#if defined ( PRISM_QUAD_TYPE )
         else if (datatype == MPI_REAL16) then

            do ipart = 1, npart
               call psmile_mg_search_2d_quad (grid_id,                       &
                           found(ipart,1)%vector, locations(ipart,1)%vector, &
                           len  (ipart,1), search_data, ipart, qtol, ierror)
               if (ierror > 0) return

               call psmile_mg_search_1d_quad (grid_id, 3,                    &
                           found(ipart,2)%vector, locations(ipart,2)%vector, &
                           search_data%search_quad(3, ipart)%vector,         &
                           search_data%dim_size(3, ipart), qtol, ierror)
               if (ierror > 0) return
            end do

#endif
         endif

      case (ndim_1d) ! 1 independent dimension
                     ! Irregular in lonlat   and vertical direction

         if (datatype == MPI_REAL) then

            call psmile_mg_search_3d_real (comp_info,      &
                        found(:, 1), locations (:, 1),     &
                        len  (:, 1), search_data, grid_id, &
                        rtol, ierror)
            if (ierror > 0) return

         else if (datatype == MPI_DOUBLE_PRECISION) then

            call psmile_mg_search_3d_dble (comp_info,      &
                        found(:, 1), locations (:, 1),     &
                        len  (:, 1), search_data, grid_id, &
                        tol, ierror)
            if (ierror > 0) return

#if defined ( PRISM_QUAD_TYPE )
         else if (datatype == MPI_REAL16) then

            call psmile_mg_search_3d_quad (comp_info,      &
                        found(:, 1), locations (:, 1),     &
                        len  (:, 1), search_data, grid_id, &
                        qtol, ierror)
            if (ierror > 0) return

#endif
         endif

      end select ! n_vec
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      end subroutine psmile_mg_search