psmile_transform_extent.F90

Go to the documentation of this file.

!------------------------------------------------------------------------
! Copyright 2010, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!------------------------------------------------------------------------

      ! the function stores the information given in shiftCount
      ! into a single integer
      ! the code is defined as follows
      ! 0000 0000 (represent the lowest 8 bits of the code)
      !       ^^^ number of shifts for the first dimension
      !      ^    sign bit for the number of shifts in the first dimension
      !  ^^^      number of shifts for the second dimension
      ! ^         sign bit for the number of shifts in the second dimension
      ! example:
      ! code = 146 ('1001 0010'B) => add 2*360 to first dimension and
      !                              subtract 1*180 from the second one
      integer function getTrCode(shiftCount)
         integer, intent(in) :: shiftCount(2)
         integer :: i

         getTrCode = 0
         do i = 2, 1, -1
            getTrCode = ishft(getTrCode, 4)
            if (shiftCount(i) < 0) getTrCode = ibset(getTrCode,3)
            getTrCode = getTrCode + abs(shiftCount(i))
         enddo
      end function getTrCode

      ! gets the shiftCount from the given code
      function getShiftCount(tr_code)
         integer, intent(in) :: tr_code
         integer :: getShiftCount(2)
         integer :: i, code

         code = tr_code
         do i = 1, 2
            getShiftCount(i) = iand(code,7)
            if (btest(code,3)) getShiftCount(i) = -1 * getShiftCount(i)
            code = ishft(code, -4)
         enddo
      end function getShiftCount

      subroutine psmile_transform_extent_cyclic (grid_type, extent, &
                              transformed, tr_codes, n_trans, ierror)
!
! !USES:
!
      use psmile_common, only : psmile_float_kind, ndim_3d, ch_id
      use psmile_grid, dummy_interface => psmile_transform_extent_cyclic

      implicit none

!
! !INPUT PARAMETERS:
!
!     prism grid type of the grid containing the extent
      integer,                  intent (in)  :: grid_type
!     Extension to be transformed and/or split
      real (psmile_float_kind), intent (in)  :: extent (2, ndim_3d)
!
! !OUTPUT PARAMETERS:
!
!     Splitted and transformed parts of the extent
      real (psmile_float_kind), intent (out) :: transformed (2, ndim_3d, 
                                                             max_num_trans_parts)
!     Code for the transformation performed
!     used by psmile_transform_extent_back to reconstruct the original range of the parts
      integer, intent (out)                  :: tr_codes (max_num_trans_parts)
!     Number of resulting parts splitted and/or transformed parts of extent
      integer, intent (out)                  :: n_trans
!     Returns the error code of psmile_Transform_extent_cyclic;
!             ierror = 0 : No error
!             ierror > 0 : Severe error (currently not used)
      integer, intent (out)                  :: ierror

!
! !LOCAL VARIABLES
!
!     control variables
      integer                  :: i, j
!     contains information on whether the grid is cyclic in the first
!     and second dimension
      logical                  :: isCyclic (2)
!     counts number of splits for a dimension
      integer                  :: splitCount
!     counts the number of shifts in each dimension that have to be performed
      integer                  :: shiftCount(2, max_num_trans_parts)
!     size of a shift for each dimension (typically [360,180])
      real (psmile_float_kind) :: shiftSize(2)
!     function declaration
      integer, external :: getTrCode
!
! !DESCRIPTION:
!
! Subroutine "psmile_Transform_extent_cyclic"
!   (*) transforms the extent of a grid of a given grid type
!       into the common coordinate system and
!   (*) splits the extent corresponding to the cyclic behaviour of
!       the given grid type
!
! The routine assumes that
! (*) longitude range is within interval (-7*180:7*180)
! (*) latitude  range is within interval (-7* 90:7* 90)
!
! and maps the data into the common_grid_range (see psmile_grid).

#ifdef VERBOSE
      print 9990, trim(ch_id), extent
      call psmile_flushstd
#endif /* VERBOSE */

      ! there is a isCyclic parameter given in the smioc,
      ! it should be taken into account as well
      ! checks whether the given grid type is cyclic in the first and/or
      ! second dimension
      isCyclic(1) = any(grids_cyclic_in_first_dim == grid_type)
      isCyclic(2) = any(grids_cyclic_in_second_dim == grid_type)

      ! computes the extent of the common grid range (typically [360,180])
      shiftSize = common_grid_range(2,1:2) - common_grid_range(1,1:2)
      shiftCount = 0
      tr_codes(:) = 0
      transformed(:,:,:) = 0
      ierror = 0

#ifdef PRISM_ASSERTION
!
!     test input data

      do i = 1, 2
         if (isCyclic(i)) then
            if ((extent(2,i) - extent(1,i)) > &
               2*(common_grid_range(2,i) - common_grid_range(1,i))) then
               print *, 'index, extent(:,index)', i, extent(:,i)
               print *, 'range is bigger than expected range: ', &
                  2*(common_grid_range(2,i) - common_grid_range(1,i))
               call psmile_assert (__FILE__, __LINE__, &
                                    'extent is too big')
            endif
            if (extent (1, i) < 7 * common_grid_range(1,i) .or. &
                extent (2, i) > 7 * common_grid_range(2,i)) then
               print *, 'index, extent', i, extent (1, i), extent (2, i)
               print *, 'range expected', 7 * common_grid_range(1,i), &
                                          7 * common_grid_range(2,i)
               call psmile_assert (__FILE__, __LINE__, &
                                   'extent is out of expected range')
            endif
            if (extent (1, i) > extent (2, i)) then
               print *, 'index, extent', i, extent (1, i), extent (2, i)
               call psmile_assert (__FILE__, __LINE__, &
                                   'upper and lower bound do not match')
            endif
         endif
      enddo

#endif

      transformed(:, :, 1) = extent(:,:)

      !
      ! 1. performs initial shifts and checks for required splits
      !
      ! for the first and second dimension
      do i = 1, 2
         if (isCyclic(i)) then
            !
            ! shift range until lower bound of extent is within the common grid range
            !
            ! shift up (+360) if necessary
            do while (transformed(1, i, 1) < common_grid_range(1,i))
               shiftCount(i,1) = shiftCount(i,1) + 1
               transformed(:,i, 1) = transformed(:,i, 1) + &
                                     shiftSize(i)
            enddo
            ! shift down (-360) if necessary
            do while (transformed(1, i, 1) >= common_grid_range(2,i))
               shiftCount(i,1) = shiftCount(i,1) - 1
               transformed(:,i, 1) = transformed(:,i, 1) - &
                                     shiftSize(i)
            enddo
         endif
      enddo

      n_trans = 1
      ! generate code for the first part
      tr_codes(1) = getTrCode(shiftCount)

      !
      ! 2. generates additional parts and their codes if required
      !
      ! for the first and second dimension
      do i = 1, 2
         if (isCyclic(i)) then

            ! initiate splitCount for current dimension
            splitCount = 0

            ! while the upper boundary exceeds the common grid range
            do while (transformed (2, i, 1 + splitCount * n_trans) > common_grid_range(2,i))

               ! the parts need to be splitted
               ! generate new parts from copies of the existing parts
               transformed(:,:,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) = &
                  transformed(:,:,1+splitCount*n_trans:(splitCount+1)*n_trans)
               shiftCount(:,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) = &
                  shiftCount(:,1+splitCount*n_trans:(splitCount+1)*n_trans)

               ! truncate upper bound of the respective dimension of the original parts
               transformed(2,i,1+splitCount*n_trans:(splitCount+1)*n_trans) = &
                  common_grid_range(2,i)

               ! shift new parts into the common range and adjust
               ! shiftCount of the new parts
               transformed(2,i,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) = &
                  transformed(2,i,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) - &
                  shiftSize(i)
               shiftCount(i,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) = &
                  shiftCount(i,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) - 1

               ! truncate lower bound of the respective dimension of new parts
               transformed(1,i,1+(splitCount+1)*n_trans:(splitCount+2)*n_trans) = &
                  common_grid_range(1,i)

               ! generate codes for new parts
               do j = 1+(splitCount+1)*n_trans, (splitCount+2)*n_trans
                  tr_codes(j) = getTrCode(shiftCount(:,j))
               enddo

               ! increase splitCount
               splitCount = splitCount + 1
            enddo

            ! add the newly created parts to the number of parts
            n_trans = n_trans * (1 + splitCount)
         endif
      enddo

#ifdef VERBOSE
      print 9980, trim(ch_id), ierror, n_trans
      call psmile_flushstd
#endif /* VERBOSE */

#ifdef VERBOSE

9990 format (1x, a, ': psmile_transform_extent_cyclic: extent', 6e13.6)
9980 format (1x, a, ': psmile_transform_extent_cyclic: eof ierror =', i3, &
                    '; n_trans =', i2)

#endif /* VERBOSE */
      end subroutine psmile_Transform_extent_cyclic
!
!--------------------------------------------------------------------------------------------------
!
      subroutine psmile_transform_extent_back (tr_codes, &
                        extents, transformed, n_trans, ierror)
!
! !USES:
!
      use psmile_common, only : psmile_float_kind, ndim_3d, ch_id
      use psmile_grid, dummy_interface => psmile_transform_extent_back

      implicit none
!
! !INPUT PARAMETERS:
!
!     Number of sub-divided extents to be transformed
      integer,                  intent (in)  :: n_trans
!     Info's to the extensions transformed
      integer,                  intent (in)  :: tr_codes (n_trans)
!     Extension to be transformed and/or subdivided
      real (psmile_float_kind), intent (in)  :: extents (2, ndim_3d, n_trans)
!
! !OUTPUT PARAMETERS:
!
!     Transformed extentions of "extents"
      real (psmile_float_kind), intent (out) :: transformed (2, ndim_3d, 
                                                             n_trans)
!     Returns the error code of psmile_Transform_extent_back;
!             ierror = 0 : No error
!             ierror > 0 : Severe error (currently not used)
      integer,                  intent (out) :: ierror
!
! !LOCAL VARIABLES
!
!     control variable
      integer :: i
!     counts the number of shifts in each dimension that have to be performed
!     (decoded from tr_code)
      integer :: shiftCount (2)
!     size of a shift for each dimension (typically [360,180])
      integer :: shiftSize(2)
!     function declaration
      interface
         function getShiftCount(tr_code)
            integer, intent(in) :: tr_code
            integer :: getShiftCount(2)
         end function getShiftCount
      end interface
!
! !DESCRIPTION:
!
! Subroutine "psmile_transform_extent_back" shifts the parts generated by
! psmile_transform_extent_cyclic back into their original coordinate range
!
#ifdef VERBOSE
      print 9990, trim(ch_id), n_trans
      call psmile_flushstd
#endif /* VERBOSE */

      ierror = 0
      transformed = extents
      ! computes the extent of the common grid range (typically [360,180])
      shiftSize = common_grid_range(2,1:2) - common_grid_range(1,1:2)

#ifdef PRISM_ASSERTION
!
!     test given input data
!
      if (any(tr_codes(:) > 255)) then
         call psmile_assert (__FILE__, __LINE__, &
                             'invalid transformation code')
      endif
#endif

      do i = 1, n_trans
         shiftCount = getShiftCount(tr_codes(i))
         transformed(1, 1:2, i) = extents(1,1:2,i) + &
                                  shiftCount * shiftSize * (-1.0)
         transformed(2, 1:2, i) = extents(2,1:2,i) + &
                                  shiftCount * shiftSize * (-1.0)
      enddo

#ifdef VERBOSE
      print 9980, trim(ch_id), ierror
      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

#ifdef VERBOSE

9990 format (1x, a, ': psmile_transform_extent_back: n_trans', i5)
9980 format (1x, a, ': psmile_transform_extent_back: eof ierror =', i3)

#endif /* VERBOSE */

      end subroutine psmile_transform_extent_back
!
!--------------------------------------------------------------------------------------------------
!
      subroutine psmile_transform_coords_db_re (tr_code_to, tr_code_from, &
                          coords_data_dble, coords_data_real, &
                          coords_size, datatype, ierror)
!
! !USES:
!
      use prism_constants
      use psmile_common, only : dble_vector, real_vector, &
                                ndim_3d, ch_id, MPI_REAL, &
                                MPI_DOUBLE_PRECISION
      use psmile_grid, dummy_interface => psmile_transform_coords_db_re

      implicit none
!
! !INPUT PARAMETERS:
!
!     Transformation code of the part in which the data is searched
      integer, intent (in) :: tr_code_to
!     Transformation code of the part which is sent and whose
!     coordinates have to be searched.
      integer, intent (in) :: tr_code_from
!     Number of coordinates to be searched
      integer, intent (in) :: coords_size (ndim_3d)
!     Datatype of input data
      integer, intent (in) :: datatype
!
! !INPUT/OUTPUT PARAMETERS:
!
!     Coordinates to be searched
      type (dble_vector), intent(inout), optional :: coords_data_dble (ndim_3d)
      type (real_vector), intent(inout), optional :: coords_data_real (ndim_3d)
!
! !OUTPUT PARAMETERS:
!
!     Returns the error code of psmile_transform_coords_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error (currently not used)
      integer, intent (out)             :: ierror
!
! !LOCAL VARIABLES
!
!     control variable
      integer :: i
!     counts the number of shifts in each dimension that have to be performed
!     (decoded from tr_code)
      integer :: shiftCount (2)
!     size of a shift for each dimension (typically [360,180])
      integer :: shiftSize(2)
!     function declaration
      interface
         function getShiftCount(tr_code)
            integer, intent(in) :: tr_code
            integer :: getShiftCount(2)
         end function getShiftCount
      end interface
!
! !DESCRIPTION:
!
! Subroutine "psmile_transform_coords_dble" transforms the
! coordinates sent by a process to be searched into the coordinate
! range of the grid in which is searched.
!
#ifdef VERBOSE
      print 9990, trim(ch_id), tr_code_to, tr_code_from
      call psmile_flushstd
#endif /* VERBOSE */

#ifdef PRISM_ASSERTION
!
!     test input data
!
      if (present(coords_data_dble) .eqv. present(coords_data_real)) then
         print *, 'coords_data_dble and coords_data_real are both', &
                  ' provided or both are missing'
         call psmile_assert (__FILE__, __LINE__, &
                              'wrong input data')
      endif
      if ((present(coords_data_dble).and.(datatype==MPI_DOUBLE_PRECISION)) .eqv. &
          (present(coords_data_real).and.(datatype==MPI_REAL))) then
         print *, 'provided data and datatype do not match'
         call psmile_assert (__FILE__, __LINE__, &
                              'wrong input data')
      endif
#endif

      ierror = 0

      if (tr_code_to /= tr_code_from) then

         ! computes the extent of the common grid range (typically [360,180])
         shiftSize = common_grid_range(2,1:2) - common_grid_range(1,1:2)

         ! gets shiftCount which is required to transform the data sent
         ! into the common grid range
         shiftCount = getShiftCount(tr_code_from)

         ! substract the shiftCount that is required to transform the
         ! local data into the common grid range
         ! this gives the shiftCount required to transform the data sent
         ! into the grid range of the local process
         shiftCount = shiftCount - getShiftCount(tr_code_to)

         ! transform the data sent
         do i = 1, 2
            select case (datatype)
               case (MPI_DOUBLE_PRECISION)
                  coords_data_dble(i)%vector(1:coords_size(i)) = &
                     coords_data_dble(i)%vector(1:coords_size(i)) + &
                     shiftCount(i) * shiftSize(i)
               case (MPI_REAL)
                  coords_data_real(i)%vector(1:coords_size(i)) = &
                     coords_data_real(i)%vector(1:coords_size(i)) + &
                     shiftCount(i) * shiftSize(i)
            end select
         enddo
      endif

#ifdef VERBOSE
      print 9980, trim(ch_id), ierror
      call psmile_flushstd
#endif /* VERBOSE */

#ifdef VERBOSE
9990 format (1x, a, ': psmile_transform_coords_db_re: ', 2i5)
9980 format (1x, a, ': psmile_transform_coords_db_re: eof ierror =', i3)
#endif /* VERBOSE */
      end subroutine psmile_transform_coords_db_re