psmile_mg_found_loc_to_3d.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_MG_found_loc_to_3d
!
! !INTERFACE:

      subroutine psmile_mg_found_loc_to_3d (search, nlev,  &
                      source_grid_type,                    &
                      found, locations, len,               &
                      virtual, virtual_cell_required,      &
                      found_3d, locations_3d, virtual_3d, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_MG_found_loc_to_3d

      implicit none
!
! !INPUT PARAMETERS:
!
      Type (Enddef_search), Intent (In)  :: search

!     Info's on coordinates to be searched
!
      Integer, Intent (In)               :: nlev

!     Number of coords to be searched in first 2 irregular directions.
!
      Integer, Intent (In)               :: source_grid_type

!     Grid type of source grid (i.e. on which the points are searched)

      Integer, Intent (In)               :: len (search%npart, *)

!     Number of coords to be searched in all directions.
!     The second dimension depends on the grid type of
!     the source grid "source_grid_type". The second dimension is:
!     1 for a source grid of type "PRISM_Irrlonlatvrt"
!           len (*)    = Number of coordinates to be searched
!     2 for a source grid of type "PRISM_Irrlonlat_Regvrt"
!           len (*, 1) = Number of coords to be searched in first 2 irregular
!                        directions.
!           len (*, 2) = Number of coords to be searched in regular z- direction.
!     3 for a source grid of type "PRISM_Reglonlatvrt"
!           len (*, i) = Number of coordinates to be searched in "i-th" direction
! |TODO
!     ! Unschoen ! Vereinheitlichen
!
      Type (integer_vector), Intent (In) :: found (search%npart, ndim_3d)
!
!     Finest level number on which a grid cell was found in the
!     method grid for point I.
!     found()%vector (i) =  1 = val_direct :
!                               Point was found on a point of the method grid.
!     found()%vector (i) = -1 = val_coupler:
!                               Point was found in a cell  of the method grid
!     Level number < -1       : Last level number on which point was found.
!     Level number = -(nlev+1): Never found (input value)
!
      Type (integer_vector), Intent (In) :: locations (search%npart, ndim_3d)
!
!     Indices of the grid cell 
!     in which the point was found in the method (corner) grid.
!
      Type (integer_vector), Intent (In) :: virtual (search%npart)

!     Code for virtual cells if point was found in virtual method cell
!     Currently available only for grids of
!     type "PRISM_Gaussreduced_regvrt".
!
      Logical,               Intent (In) :: virtual_cell_required
!
!     Is the virtual cell info available and is it required to 
!     generate the 3d virtual cell info ?
!
! !OUTPUT PARAMETERS:
!
      Type (integer_vector), Intent(Out) :: found_3d     (search%npart)

!     Finest level number on which a grid cell was found in the
!     method grid for point I stored as 3d arrays.

      Type (integer_vector), Intent(Out) :: locations_3d (search%npart)
!
!     Indices of the grid cell 
!     in which the point was found in the method (corner) grid stored
!     as 3d arrays.
!
      Type (integer_vector), Intent(Out) :: virtual_3d (search%npart)
!
!     Code for the virtual cell
!     in which the point was found in the method (corner) grid stored
!     as 3d arrays. Only set if "virtual_cell_required" == .true.
!
      Integer, Intent (Out)              :: ierror

!     Returns the error code of PSMILe_MG_found_loc_to_3d;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
! indl = Index of LonLat values in "found, locations"
! indz = Index of Vert   values in "found, locations"
!
      Integer, Parameter           :: indl = 1
      Integer, Parameter           :: indz = 2
!
      Integer                      :: i, j, k, ipart, npart
      Integer                      :: search_grid_type
      Integer                      :: not_found
!
!     ... for 
!
      Integer                      :: len2, len3, nprev
      Integer                      :: leni, lenj, lenk
      Integer                      :: leni_mask, lenj_mask, lenk_mask
      Integer                      :: ii, jj
!
!     ... for error parameters
!
      Integer, parameter           :: nerrp = 2
      Integer                      :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_MG_found_loc_to_3d" transforms the input fields
! "found" and "locations" (containing the results of the MG search)
! into the fully 3d versions of that fields.
! Additionally, the mask "search_mask" of the points to be searched
! is applied if required.
!
! TODO: irrlonlat_regvrt and gauss_reduced should be handled
!       identically
!
! !REVISION HISTORY:
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.05    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_mg_found_loc_to_3d.F90 2788 2010-11-30 14:34:07Z hanke $
! $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
     '$Id: psmile_mg_found_loc_to_3d.F90 2788 2010-11-30 14:34:07Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */

      npart = search%npart
      search_grid_type = search%grid_type
!
!===>  Allocate full 3d vectors of "found" and "location"
!      locations_3d koennte man sich sparen, wenn man die locations fuer
!      jeden Index einzeln als vector abspeichert.
!
         do ipart = 1, npart

         len3 = (search%range(2,1,ipart)-search%range(1,1,ipart)+1) * &
                (search%range(2,2,ipart)-search%range(1,2,ipart)+1) * &
                (search%range(2,3,ipart)-search%range(1,3,ipart)+1)

#ifdef PRISM_ASSERTION
         if (source_grid_type == PRISM_Reglonlatvrt) then
            if (search_grid_type == PRISM_Irrlonlatvrt) then
                  do i = 1, ndim_3d
                  if (len3 /= len (ipart,i)) then
                     call psmile_assert (__FILE__, __LINE__, &
                                         "len3 /= len(ipart,i)")
                  endif 
                  end do
            else if (search_grid_type == PRISM_Irrlonlat_Regvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,3)) then
                     print *, 'len3, len', len3, len (ipart,1:ndim_3d)
                     call psmile_assert (__FILE__, __LINE__, &
                                    "len3 /= len(ipart,1)*len(ipart,3)")
                  endif 
            else if (search_grid_type == PRISM_Reglonlatvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,2)*len(ipart,3)) then
                     print *, 'len3, len', len3, len (ipart,1:ndim_3d)
                     call psmile_assert (__FILE__, __LINE__, &
                                 "len3 /= Product(len(ipart,1:ndim3d))")
                  endif 
            else if (search_grid_type == PRISM_Gaussreduced_regvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,3)) then
                     print *, 'len3, len', len3, len (ipart,1:ndim_3d)
                     call psmile_assert (__FILE__, __LINE__, &
                                    "len3 /= len(ipart,1)*len(ipart,3)")
                  endif 
           endif
!
         else if (source_grid_type == PRISM_Irrlonlat_Regvrt) then
            if (search_grid_type == PRISM_Irrlonlatvrt) then
                  do i = 1, 2
                  if (len3 /= len (ipart,i)) then
                     call psmile_assert (__FILE__, __LINE__, &
                                         "len3 /= len(ipart,i)")
                  endif 
                  end do
            else if (search_grid_type == PRISM_Gaussreduced_regvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,2)) then
                     print *, 'len3, len', len3, len (ipart,1:2)
                     call psmile_assert (__FILE__, __LINE__, &
                                  "len3 /= len(ipart,1)*len(ipart,2)")
                  endif 
            else ! if (search_grid_type == PRISM_Irrlonlat_Regvrt .or.
                 !     search_grid_type == PRISM_Reglonlatvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,2)) then
                     print *, 'len3, len', len3, len (ipart,1:2)
                     call psmile_assert (__FILE__, __LINE__, &
                                  "len3 /= len(ipart,1)*len(ipart,2)")
                  endif 
            endif
!
         else if (source_grid_type == PRISM_Gaussreduced_Regvrt) then
            if (search_grid_type == PRISM_Irrlonlatvrt) then
                  do i = 1, 2
                  if (len3 /= len (ipart,i)) then
                     call psmile_assert (__FILE__, __LINE__, &
                                                "len3 /= len(ipart,i)")
                  endif 
                  end do
            else if (search_grid_type == PRISM_Gaussreduced_regvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,2)) then
                     print *, 'len3, len', len3, len (ipart,1:2)
                     call psmile_assert (__FILE__, __LINE__, &
                                  "len3 /= len(ipart,1)*len(ipart,3)")
                  endif 
            else ! if (search_grid_type == PRISM_Irrlonlat_Regvrt .or.
                 !     search_grid_type == PRISM_Reglonlatvrt) then
                  if (len3 /= len (ipart,1)*len(ipart,2) ) then
                     call psmile_assert (__FILE__, __LINE__, &
                                   "len3 /= len(ipart,1)*len(ipart,2)")
                  endif
            endif
!
         else if (source_grid_type == PRISM_Irrlonlatvrt) then
            if (len3 /= len (ipart,1)) then
               call psmile_assert (__FILE__, __LINE__, &
                                   "len3 /= len(ipart,1)")
            endif 
         endif
!
         if (virtual_cell_required .and. &
             source_grid_type /= PRISM_Gaussreduced_regvrt) then

            call psmile_assert (__FILE__, __LINE__, &
               "virtual_cell_required only available for PRISM_Gaussreduced_regvrt")
         endif
#endif

         Allocate (found_3d(ipart)%vector(len3), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = len3
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'found_3d(ipart)%vector', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
         Allocate (locations_3d(ipart)%vector(ndim_3d*len3), STAT = ierror)
         if ( ierror > 0 ) then
            ierrp (1) = ierror
            ierrp (2) = len3 * ndim_3d
            ierror = PRISM_Error_Alloc
            call psmile_error ( ierror, 'locations_3d(ipart)%vector', &
                                ierrp, 2, __FILE__, __LINE__ )
            return
         endif
!
         if (virtual_cell_required) then
            Allocate (virtual_3d(ipart)%vector(len3), STAT = ierror)
            if ( ierror > 0 ) then
               ierrp (1) = ierror
               ierrp (2) = len3
               ierror = PRISM_Error_Alloc
               call psmile_error ( ierror, 'virtual(ipart)%vector', &
                                   ierrp, 2, __FILE__, __LINE__ )
               return
            endif
         endif
         end do ! ipart
!
!-----------------------------------------------------------------------
!      Source grid is regular in all directions;
!      i.e. all target points are searched for all directions
!      separately.
!-----------------------------------------------------------------------
!
      if (source_grid_type == PRISM_Reglonlatvrt) then
!
         if (search_grid_type == PRISM_Irrlonlatvrt) then
!
!===> ... Target grid is of type PRISM_Irrlonlatvrt;
!         i.e. target grid is fully irregular in all directions and
!              points to be searched are stored as 3d array.
!         len(ipart,1) == len(ipart,2) == len(ipart,3)
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
                  found_3d(ipart)%vector(i) = &
                     min (found(ipart,1)%vector(i), &
                          found(ipart,2)%vector(i), &
                          found(ipart,3)%vector(i))
                  end do ! i
               end do ! ipart
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+1) = &
                     locations(ipart,1)%vector(i)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+2) = &
                     locations(ipart,2)%vector(i)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+3) = &
                     locations(ipart,3)%vector(i)
                  end do ! i

               end do ! ipart

         else if (search_grid_type == PRISM_Irrlonlat_Regvrt) then
!
!===> ... Target grid is of type PRISM_Irrlonlat_Regvrt;
!         i.e. target grid is irregular in LonLat directions and
!              regular in z-direction.
!              Points to be searched are stored as 2d array and 1d
!              array.
!         len(ipart,1) = Number of points in Lon and Lat direction (2d)
!         len(ipart,2) = Number of points in Lon and Lat direction (2d)
!         len(ipart,3) = Number of points in z-direction (1d)
!
               do ipart = 1, npart
                  do k = 1, len(ipart,3)
                  nprev = (k-1)*len(ipart,1)
!cdir vector
                     do i = 1, len(ipart,1)
                     found_3d(ipart)%vector(nprev+i) = &
                        min (found(ipart,1)%vector(i), &
                             found(ipart,2)%vector(i), &
                             found(ipart,3)%vector(k))
                     end do ! i
                  end do ! k
               end do ! ipart
!
               do ipart = 1, npart
                  do k = 1, len(ipart,3)
                  nprev = (k-1)*len(ipart,1) * ndim_3d
!cdir vector
                     do i = 1, len(ipart,1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,1)%vector(i)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                        locations(ipart,2)%vector(i)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                        locations(ipart,3)%vector(k)
                     end do ! i

                  end do ! k
               end do ! ipart
!
         else if (search_grid_type == PRISM_Gaussreduced_Regvrt) then
!
!===> ... Target grid is of type PRISM_Gaussreduced_Regvrt;
!         i.e. target grid is irregular in LonLat directions and
!              regular in z-direction.
!              Points to be searched are stored as 2d array and 1d
!              array.
!         len(ipart,1) = Number of points in Lon and Lat direction (2d)
!         len(ipart,2) = Number of points in Lon and Lat direction (2d)
!         len(ipart,3) = Number of points in z-direction (1d)
!
               do ipart = 1, npart

                  do k = 1, len(ipart,3)
                  nprev = (k-1)*len(ipart,1)
!cdir vector
                     do i = 1, len(ipart,1)
                     found_3d(ipart)%vector(nprev+i) = &
                        min (found(ipart,1)%vector(i), &
                             found(ipart,2)%vector(i), &
                             found(ipart,3)%vector(k))
                     end do ! i
                  end do ! k
               end do ! ipart
!
               do ipart = 1, npart
                  do k = 1, len(ipart,3)
                  nprev = (k-1)*len(ipart,1) * ndim_3d
!cdir vector
                     do i = 1, len(ipart,1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,1)%vector(i)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                        locations(ipart,2)%vector(i)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                        locations(ipart,3)%vector(k)
                     end do ! i

                  end do ! k
               end do ! ipart
!
         else if (search_grid_type == PRISM_Reglonlatvrt) then
!
!===> ... Target grid is of type PRISM_Reglonlatvrt;
!         i.e. target grid is regular in all directions.
!              Points to be searched are stored as 1d arrays.
!
! ? Loops splitten und konstante Teile einzeln zuweisen ?
!
               do ipart = 1, npart
               len2 = len(ipart,1) * len(ipart,2)
!
                  do k = 1, len(ipart,3)
                     do j = 1, len(ipart,2)
                     nprev = (k-1)*len2 + (j-1)*len(ipart,1)
!                    min1 = min (found(ipart,2)%vector(j), &
!                                found(ipart,3)%vector(k))
!cdir vector
                        do i = 1, len(ipart,1)
                        found_3d(ipart)%vector(nprev+i) = &
                           min (found(ipart,1)%vector(i), &
                                found(ipart,2)%vector(j), &
                                found(ipart,3)%vector(k))
                        end do ! i
                     end do ! j
                  end do ! k
               end do ! ipart
!
               do ipart = 1, npart
               len2 = len(ipart,1) * len(ipart,2)
!
                  do k = 1, len(ipart,3)
                     do j = 1, len(ipart,2)
                     nprev = ((k-1)*len2 + (j-1)*len(ipart,1)) * ndim_3d
!cdir vector
                        do i = 1, len(ipart,1)
                        locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                           locations(ipart,1)%vector(i)
                        locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                           locations(ipart,2)%vector(j)
                        locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                           locations(ipart,3)%vector(k)
                        end do ! i

                     end do ! j
                  end do ! k
               end do ! ipart
         else
!
!===> ... Unknown target grid type
!
            ierrp (1) = search_grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'Unknown search grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
            return
         endif
!
!-----------------------------------------------------------------------
!      Source grid is irregular in LonLat direction and
!      regular in z direction.
!      i.e. all target points are searched for lonlat direction
!      and z direction separately.
!-----------------------------------------------------------------------
!
      else if ( source_grid_type == PRISM_Irrlonlat_Regvrt ) then
!
         if (search_grid_type == PRISM_Irrlonlatvrt) then
!
!===> ... Target grid is of type PRISM_Irrlonlatvrt;
!         i.e. target grid is fully irregular in all directions and
!              points to be searched are stored as 3d array.
!         len(ipart,1) == len(ipart,2)
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
                  found_3d(ipart)%vector(i) = &
                     min (found(ipart,indl)%vector(i), &
                          found(ipart,indz)%vector(i))
                  end do ! i
               end do ! ipart
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+1) = &
                     locations(ipart,indl)%vector((i-1)*ndim_2d+1)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+2) = &
                     locations(ipart,indl)%vector((i-1)*ndim_2d+2)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+3) = &
                     locations(ipart,indz)%vector(i)
                  end do ! i

               end do ! ipart

         else if (search_grid_type == PRISM_Irrlonlat_Regvrt .or. &
                  search_grid_type == PRISM_Reglonlatvrt) then
!
!===> ... Target grid is of type PRISM_Irrlonlat_Regvrt;
!         i.e. target grid is irregular in LonLat directions and
!              regular in z-direction.
!              Points to be searched are stored as 2d array and 1d
!              array.
!         or Target grid is of type PRISM_Reglonlatvrt;
!         i.e. target grid is regular in all directions.
!              Points to be searched were generated as 2d array and
!              1d array.
!         len(ipart,1) = Number of points in Lon and Lat direction (2d)
!         len(ipart,2) = Number of points in z-direction (1d)
!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1)
!cdir vector
                     do i = 1, len(ipart,1)
                     found_3d(ipart)%vector(nprev+i) =    &
                        min (found(ipart,indl)%vector(i), &
                             found(ipart,indz)%vector(k))
                     end do ! i
                  end do ! k
               end do ! ipart

!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1) * ndim_3d
!cdir vector
                     do i = 1, len(ipart,1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+2)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                        locations(ipart,indz)%vector(k)
                     end do ! i

                  end do ! k
               end do ! ipart

         else if (search_grid_type == PRISM_Gaussreduced_Regvrt) then
!
!===> ... Target grid is of type PRISM_Gaussreduced_Regvrt;
!         i.e. target grid is irregular in LonLat directions and
!              regular in z-direction.
!              Points to be searched are stored as 2d array and 1d
!              array.
!         len(ipart,1) = Number of points in Lon and Lat direction (2d)
!         len(ipart,2) = Number of points in Lon and Lat direction (2d)
!         len(ipart,3) = Number of points in z-direction (1d)
!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1)
!cdir vector
                     do i = 1, len(ipart,1)
                     if ( found(ipart,indl)%vector(i) > 1 .or. &
                          found(ipart,indz)%vector(k) > 1 ) then
                        found_3d(ipart)%vector(nprev+i) = -(nlev+1)
                     else
                        found_3d(ipart)%vector(nprev+i) = &
                           min (found(ipart,indl)%vector(i), &
                                found(ipart,indz)%vector(k))
                     endif
                     end do ! i
                  end do ! k
               end do ! ipart
!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1) * ndim_3d
!cdir vector
                     do i = 1, len(ipart,1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+2)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                        locations(ipart,indz)%vector(k)
                     end do ! i

                  end do ! k
               end do ! ipart

         else
!
!===> ... Unknown target grid type
!
            ierrp (1) = search_grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'Unknown search grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
            return
         endif
!
!-----------------------------------------------------------------------
!     Source grid is GaussReduced grid
!-----------------------------------------------------------------------
!
      else if (source_grid_type == PRISM_Gaussreduced_Regvrt ) then
!
         if (search_grid_type == PRISM_Irrlonlatvrt) then
!
!===> ... Target grid is of type PRISM_Irrlonlatvrt;
!         i.e. target grid is fully irregular in all directions and
!              points to be searched are stored as 3d array.
!         len(ipart,1) == len(ipart,2)
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
                  found_3d(ipart)%vector(i) = &
                     min (found(ipart,indl)%vector(i), &
                          found(ipart,indz)%vector(i))
                  end do ! i
               end do ! ipart
!
               do ipart = 1, npart
!cdir vector
                  do i = 1, len(ipart,1)
#ifdef STORE_LOCATION_AS_2D
                  locations_3d(ipart)%vector((i-1)*ndim_3d+1) = &
                     locations(ipart,indl)%vector((i-1)*ndim_2d+1)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+2) = &
                     locations(ipart,indl)%vector((i-1)*ndim_2d+2)
#else
                  locations_3d(ipart)%vector((i-1)*ndim_3d+1) = &
                     locations(ipart,indl)%vector(i)
                  locations_3d(ipart)%vector((i-1)*ndim_3d+2) = 1
#endif
                  locations_3d(ipart)%vector((i-1)*ndim_3d+3) = &
                     locations(ipart,indz)%vector(i)
                  end do ! i

               end do ! ipart
!
               if (virtual_cell_required) then
!                 Currently virtual cell info only available for LonLat
!                 part of GaussReduce Grid
!
                  do ipart = 1, npart
                  virtual_3d(ipart)%vector(:) = virtual(ipart)%vector(:)
                  end do ! ipart
               endif

         else if (search_grid_type == PRISM_Irrlonlat_Regvrt .or. &
                  search_grid_type == PRISM_Reglonlatvrt     .or. &
                  search_grid_type == PRISM_Gaussreduced_Regvrt ) then
!
!===> ... Target grid is of type PRISM_Irrlonlat_Regvrt;
!         i.e. target grid is irregular in LonLat directions and
!              regular in z-direction.
!              Points to be searched are stored as 2d array and 1d
!              array.
!         or Target grid is of type PRISM_Reglonlatvrt;
!         i.e. target grid is regular in all directions.
!              Points to be searched were generated as 2d array and
!              1d array.
!         len(ipart,1) = Number of points in Lon and Lat direction (2d)
!         len(ipart,2) = Number of points in z-direction (1d)
!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1)
!cdir vector
                     do i = 1, len(ipart,1)
                     found_3d(ipart)%vector(nprev+i) = &
                        min (found(ipart,indl)%vector(i), &
                             found(ipart,indz)%vector(k))
                     end do ! i
                  end do ! k
               end do ! ipart
!
               do ipart = 1, npart
                  do k = 1, len(ipart,2)
                  nprev = (k-1)*len(ipart,1) * ndim_3d
!cdir vector
                     do i = 1, len(ipart,1)
#ifdef STORE_LOCATION_AS_2D
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+1)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = &
                        locations(ipart,indl)%vector( (i-1)*ndim_2d+2)
#else
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+1) = &
                        locations(ipart,indl)%vector(i)
                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+2) = 1
#endif

                     locations_3d(ipart)%vector(nprev+(i-1)*ndim_3d+3) = &
                        locations(ipart,indz)%vector(k)
                     end do ! i

                  end do ! k
               end do ! ipart
!
               if (virtual_cell_required) then
!                 Currently virtual cell info only available for LonLat
!                 part of GaussReduce Grid
!
                  do ipart = 1, npart
                     do k = 1, len(ipart,2)
                     nprev = (k-1)*len(ipart,1)
                     virtual_3d(ipart)%vector(nprev+1:nprev+len(ipart,1)) = &
                        virtual(ipart)%vector(      1:      len(ipart,1))
                     end do ! k
                  end do ! ipart
               endif
!
         else
!
!===> ... Unknown target grid type
!
            ierrp (1) = search_grid_type
            ierror = PRISM_Error_Internal

            call psmile_error ( ierror, 'Unknown search grid type', &
                                ierrp, 1, __FILE__, __LINE__ )
            return
         endif
!
!-----------------------------------------------------------------------
!      Source grid is irregular in all directions
!      i.e. all target points are searched in all 3 directions.
!      Routine is not necessary
!-----------------------------------------------------------------------
!
      else if (source_grid_type == PRISM_Irrlonlatvrt) then
!
         print 9970, trim(ch_id)
!
            do ipart = 1, npart
            found_3d(ipart)%vector(:) = found(ipart,1)%vector(:)
            end do
!
            do ipart = 1, npart
            locations_3d(ipart)%vector(:) = locations(ipart,1)%vector(:)
            end do
!
!-----------------------------------------------------------------------
!     Unknown source grid type
!-----------------------------------------------------------------------
!
      else
         ierrp (1) = source_grid_type
         ierror = PRISM_Error_Internal

         call psmile_error ( ierror, 'Unknown source grid type', &
                             ierrp, 1, __FILE__, __LINE__ )
         return

      endif
!
!-----------------------------------------------------------------------
!     Apply mask on points to be searched if required;
!     i.e. set "found_3d" to -(nlev+1) ("not found")
!          if search mask is not true.
!-----------------------------------------------------------------------
!
      if (Associated(search%search_mask)) then

         not_found = - (nlev+1)
         do ipart = 1, npart

            if ( search%msg_intersections%requires_conserv_remap /= PSMILe_conserv2D .and. &
                 search%msg_intersections%requires_conserv_remap /= PSMILe_conserv3D ) then

               len3 = (search%range(2,1,ipart)-search%range(1,1,ipart)+1) * &
                      (search%range(2,2,ipart)-search%range(1,2,ipart)+1) * &
                      (search%range(2,3,ipart)-search%range(1,3,ipart)+1)
!cdir vector
               do i = 1, len3
                  if (.not. search%search_mask(ipart)%vector(i)) then
                     found_3d(ipart)%vector(i) = not_found
                  endif
               end do

            else

               leni =  search%range(2,1,ipart)-search%range(1,1,ipart)+1
               lenj =  search%range(2,2,ipart)-search%range(1,2,ipart)+1
               lenk =  search%range(2,3,ipart)-search%range(1,3,ipart)+1

               if ( search_grid_type == PRISM_Gaussreduced_Regvrt ) then

                  if ( search%msg_intersections%requires_conserv_remap == PSMILe_conserv2D ) then
                     leni_mask = (search%range(2,1,ipart)-search%range(1,1,ipart))/2 + 1
                     lenj_mask = 1
                     lenk_mask = search%range(2,3,ipart)-search%range(1,3,ipart)+1
                  else
                     leni_mask = (search%range(2,1,ipart)-search%range(1,1,ipart))/2 + 1
                     lenj_mask = 1
                     lenk_mask = search%range(2,3,ipart)-search%range(1,3,ipart)
                  endif
 
                  do k = 1, lenk_mask
                     do j = 1, lenj_mask
                        do i = 1, leni_mask
                           ii = (k-1)*lenj_mask*leni_mask+(j-1)*leni_mask+i
                           if (.not. search%search_mask(ipart)%vector(ii)) then
                              jj = (k-1)*lenj*leni+(j-1)*leni+i
                              found_3d(ipart)%vector(jj) = not_found
                           endif
                        enddo
                     enddo
                  enddo

               else

                  if ( search%msg_intersections%requires_conserv_remap == PSMILe_conserv2D ) then
                     leni_mask = search%range(2,1,ipart)-search%range(1,1,ipart)
                     lenj_mask = search%range(2,2,ipart)-search%range(1,2,ipart)
                     lenk_mask = search%range(2,3,ipart)-search%range(1,3,ipart)+1
                  else
                     leni_mask = search%range(2,1,ipart)-search%range(1,1,ipart)
                     lenj_mask = search%range(2,2,ipart)-search%range(1,2,ipart)
                     lenk_mask = search%range(2,3,ipart)-search%range(1,3,ipart)
                  endif

                  do k = 1, lenk_mask
                     do j = 1, lenj_mask
                        do i = 1, leni_mask
                           ii = (k-1)*lenj_mask*leni_mask+(j-1)*leni_mask+i
                           if (.not. search%search_mask(ipart)%vector(ii)) then
                              jj = (k-1)*lenj*leni+(j-1)*leni+i
                              found_3d(ipart)%vector(jj) = not_found
                           endif
                        enddo
                     enddo
                  enddo

               endif

            endif

         end do
      endif
!
!-----------------------------------------------------------------------
!     All done
!-----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990  format (1x, a, ': psmile_mg_found_loc_to_3d:')
9980  format (1x, a, ': psmile_mg_found_loc_to_3d: eof ierror =', i4)
!
9970  format (/1x, a, ': ### WARNING PSMILe_MG_found_loc_to_3d called ', &
                      'for 3d arrays !!!',                               &
              /1x)

      end subroutine PSMILe_MG_found_loc_to_3d