psmile_transform_gauss2.F90

Go to the documentation of this file.

       subroutine psmile_transform_gauss2 ( search, glen, map_shape,  &
                           map, locations, found, glocations, gfound, &
                           cell_based_switch, nlev, grid_valid_shape, &
                           ierror )
!
! !USES:
!
      use PRISM_constants
      use PSMILe, dummy_interface => PSMILe_Transform_gauss2
#ifdef DEBUG_TRACE
  use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Type (Enddef_search)                :: search

!     Info's on coordinates to be searched

      Integer, Intent (In)                :: glen (search%npart, 2)

!     Number of coords to be searched in all directions.
!     glen(*,1) = Number in lonlat/GaussReduced direction
!     glen(*,2) = Number in vertical direction

      Integer, Intent (In)                :: map_shape(ndim_3d)
      Integer, Intent (In)                :: map (map_shape(1), 
                                                  map_shape(2), 
                                                  map_shape(3))
!
!     Mapping from regular indices in lon and lat direction of regular grid
!     into (sequential) lonlat directions of Gaussian Reduced Grid
!     TODO: Should become indices in lon direction for current longitude
!
      Type (integer_vector), Intent (In)  :: locations (search%npart, ndim_3d)
!
!     Input:  Indices of the grid cell in which the point was found
!             in the auxillary grid.
!             Note: The auxiliary grid is of type PRISM_Reglonlatvrt.

      Type (integer_vector), Intent (In)  :: found (search%npart, ndim_3d)
!
!     Input:  Finest level number on which a grid cell was found
!             in the auxillary grid for point I mapped onto the gauss grid.
!             Level number = -(nlev+1): Never found
!
      Logical, Intent (In)                :: cell_based_switch
!
!     Indication that we have changed from cell-based to point-based search 
!
      Integer, Intent (In)                :: nlev
!
!     Maximum number of levels in the multigrid hierarchy. Needed to set the
!     gfound array to "not found" for points that were found in the auxiliary
!     grid but not on the actual reduced gauss grid.
!
      Integer, Intent (In)                :: grid_valid_shape(2)
!
!     Defines the valid coordinate range for the actual reduced gauss grid.
!
! !INPUT/OUTPUT PARAMETERS:
!
!
! !OUTPUT PARAMETERS:

      Type (integer_vector), Intent (InOut) :: glocations (search%npart, 2)
!
!     Indices of the grid cell of the point that was found
!     in the auxillary grid mapped onto the gauss grid.
!     The latitude  values are in the range of 1:nbr_lats
!     The longitude values are in the range of 1:points_per_lat(lat)
!     glocations (*,1)%vector is a 1-dimensional vector in GaussReduced
!     grid notation.

      Type (integer_vector), Intent (InOut) :: gfound (search%npart, 2)
!
!     Finest level number on which a grid cell was found
!     in the auxillary grid for point I.
!     Level number = -(nlev+1): Never found
!
      Integer,               Intent (Out) :: ierror

!     Returns the error code of PSMILe_Transform_gauss2_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
! indl = Index of LonLat/GaussReduced values in "gfound, glocations"
! indz = Index of Vert   values in "gfound, glocations"
!
      Integer, Parameter           :: indl = 1
      Integer, Parameter           :: indz = 2

      Integer                      :: i, ipart, j, npart, n
      Integer                      :: ii, jlat
!
!     ... for error parameters
!
      Integer, parameter           :: nerrp = 2
      Integer                      :: ierrp (nerrp)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Transform_gauss2" transforms the locations
! and found information from the auxillary grid onto the Gauss
! grid.
!
! !REVISION HISTORY:
!   Date      Programmer   Description
! ----------  ----------   -----------
! 07.01.05    R. Redler    created
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_transform_gauss2.F90 2966 2011-02-18 09:47:30Z hanke $
! $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_transform_gauss2.F90 2966 2011-02-18 09:47:30Z hanke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */
!
      npart = search%npart

!===> Set z-coordinate (independent on supported grid types)

      do ipart = 1, npart
         gfound (ipart,indz)%vector     = found(ipart,3)%vector
         glocations (ipart,indz)%vector = locations(ipart,3)%vector
      enddo
!
! Transform locations and found array depending on the target grid type
!
      select case ( search%grid_type )

      case ( PRISM_irrlonlat_regvrt)
!
!===> ... Generate locations in gaussian reduced grid
!
! Indication that we have changed from cell-based to point-based search
! old       if ( size( locations(ipart,1)%vector) /= &
!                size(glocations(ipart,1)%vector) )

         do ipart = 1, search%npart
!
!===> ... Generate found vector in lonlat directon
!         out of both partial found vectors in lon and lat direction
!
            gfound (ipart, 1)%vector(:) = &
               min (found (ipart, 1)%vector(:), found (ipart, 2)%vector(:))

            ii = 1
            do j = search%range(1,2,ipart), search%range(2,2,ipart)
!cdir vector
               do i = search%range(1,1,ipart), search%range(2,1,ipart)

                  n = psmile_transform_index_3d_to_1d ((/i, j, 1/), &
                                                       search%range(:,:,ipart))

                  if ( gfound(ipart,1)%vector(ii) == 1 ) then
                     glocations(ipart,1)%vector(ii) = &
                        map(locations(ipart,1)%vector(n),locations(ipart,2)%vector(n), 1)
#ifdef DEBUG_TRACE
                     if ( i == ictl_ind(1) .and. j == ictl_ind(2)) then
                        print 8990, trim(ch_id), locations(ipart,1)%vector(n), &
                                                 locations(ipart,2)%vector(n)
                     endif
8990 format (1x, a, ': psmile_transform_gauss2 ictl location on aux grid: ', 2i10)
#endif
                  else
                     glocations(ipart,1)%vector(ii) = PSMILe_undef
                  endif
                  ii = ii + 1
               end do ! i
            end do ! j
         end do !  ipart

      case ( PRISM_reglonlatvrt)
!
!===> ... Generate found vector in lonlat directon
!         out of both partial found vectors in lon and lat direction
!         The partial vectors are 1d vectors.
!         The generated vector is a 1d vector.
!
         do ipart = 1, search%npart
            ii = 1
            do jlat = 1, search%range(2,2,ipart)-search%range(1,2,ipart)+1
!cdir vector
               do i = 1, search%range(2,1,ipart)-search%range(1,1,ipart)+1

                  gfound(ipart,1)%vector(ii) = &
                     min(found(ipart,1)%vector(i), found(ipart,2)%vector(jlat))

                  if ( gfound(ipart,1)%vector(ii) == 1 ) then
                     glocations(ipart,1)%vector(ii) =        &
                          map(locations(ipart,1)%vector(i),  &
                              locations(ipart,2)%vector(jlat), 1)
                  else
                     glocations(ipart,1)%vector(ii) = PSMILe_undef
                  endif
                  ii = ii + 1
               end do ! i
            end do ! jlat
         end do ! ipart

      case ( PRISM_Gaussreduced_regvrt, PRISM_Irrlonlatvrt)

         do ipart = 1, search%npart
            gfound (ipart, 1)%vector(:) = &
               min (found (ipart, 1)%vector(:), found (ipart, 2)%vector(:))

!cdir vector
            do i = 1, glen(ipart,1)
               if ( gfound(ipart,1)%vector(i) == 1 ) then
                  glocations(ipart,1)%vector(i) = &
                     map(locations(ipart,1)%vector(i),locations(ipart,2)%vector(i), 1)
               else
                  glocations(ipart,1)%vector(i) = PSMILe_undef
               endif
            enddo ! i
         enddo ! ipart

!     case ( PRISM_Irrlonlatvrt)
!     wuerde einen vollen 3d-vector brauchen
!     wird in einer anderen Routine gemacht ?!/

      case default

         ierrp (1) = search%grid_type
         ierror = PRISM_Error_Internal
         call psmile_error ( ierror, 'unsupported search grid type', &
              ierrp, 1, __FILE__, __LINE__ )
         return

      end select

! there may be points, that were found on the auxiliary grid but are not on the actual
! reduced gauss grid. For these points, the gfound vector needs to be adjusted accordingly.

      do ipart = 1, search%npart
        where (glocations(ipart,1)%vector(:) < grid_valid_shape(1) .or. &
                glocations(ipart,1)%vector(:) > grid_valid_shape(2))

               gfound(ipart,1)%vector(:) = -(nlev+1)
         end where
      enddo

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

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

      end subroutine psmile_transform_gauss2