psmile_bbcells_pole_dble.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2010, DKRZ, Hamburg, Germany.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: psmile_bbcells_pole_dble
!
! !INTERFACE:

      subroutine psmile_bbcells_pole_dble (              &
                     coords_shape, coords_x, coords_y,   &
                     corner_shape, sub_range,            &
                     chmin_x, chmax_x, chmin_y, chmax_y, &
                     midp_x, midp_y, pole_array, period, &
                     ierror)
!
! !USES:
!
      use psmile, dummy_interface => psmile_bbcells_pole_dble
      use psmile_reallocate
      use psmile_common, only : ndim_2d, ndim_3d
      use psmile_grid, only : psmile_transform_cell_cyclic
      use prism, only : prism_undefined

      implicit none
!
! !INPUT PARAMETERS:
!
      integer,           intent (in)  :: coords_shape (2, ndim_2d)
!     Dimension of (method/grid) coordinate "coords"
      double precision,  intent (in)  :: 
         coords_x(coords_shape(1,1):coords_shape(2,1),  
                  coords_shape(1,2):coords_shape(2,2)), 
         coords_y(coords_shape(1,1):coords_shape(2,1),  
                  coords_shape(1,2):coords_shape(2,2))
!     Coordinates of the method
      integer,           intent (in)    :: corner_shape (2, ndim_3d)
!     Dimension of the array "corners" containing the corner points
!     of the control volume (for a single coordinate).
      integer,           intent (in)    :: sub_range (2, ndim_2d)
!     Definintion of the subgrid for which the (bounding) cells should
!     be computed.
      integer,           intent (in)    :: pole_array(:)
!     List of corner cells covering the poles
      double precision, intent(in)      :: period
!     period of grid coordinate range (typically 360 degree)

!
! !INOUTPUT PARAMETERS:
!
      double precision,  intent (inout) :: chmin_x (sub_range(1,1)-1:sub_range(2,1), 
                                                    sub_range(1,2)-1:sub_range(2,2))
!     Minimum value of bounding box of method cell
      double precision,  intent (inout) :: chmax_x (sub_range(1,1)-1:sub_range(2,1), 
                                                    sub_range(1,2)-1:sub_range(2,2))
!     Maximum value of bounding box of method cell
      double precision,  intent (inout) :: chmin_y (sub_range(1,1)-1:sub_range(2,1), 
                                                    sub_range(1,2)-1:sub_range(2,2))
!     Minimum value of bounding box of method cell
      double precision,  intent (inout) :: chmax_y (sub_range(1,1)-1:sub_range(2,1), 
                                                    sub_range(1,2)-1:sub_range(2,2))
!     Maximum value of bounding box of method cell
      double precision,  intent (inout) :: midp_x(sub_range(1,1)-1:sub_range(2,1), 
                                                  sub_range(1,2)-1:sub_range(2,2))
!     Mid point of bounding box of method cell
      double precision,  intent (inout) :: midp_y(sub_range(1,1)-1:sub_range(2,1), 
                                                  sub_range(1,2)-1:sub_range(2,2))
!     Mid point of bounding box of method cell

!
! !OUTPUT PARAMETERS:
!
      integer,           intent (out)   :: ierror

!     Returns the error code of psmile_bbcells_pole_dble
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
! Loop variable
      integer                      :: n, i, j
! Index for accessing arrays
      integer                      :: ijk(3)
! For computation of a new midpoint
      type(point_dble)             :: cell(4), temp_cell(4), midpoint
! Method grid pole array
      integer, pointer             :: method_pole_array_x(:), method_pole_array_y(:)
      integer                      :: method_pole_array_size
      integer                      :: method_pole_array_allocated_size
! Error parameters
      integer, parameter           :: nerrp = 2
      integer                      :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "psmile_bbcells_pole_dble" extents the bounding boxes
! of cells covering the north or south pole to 90 or -90
! respectively
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 14.09.10    M. Hanke     created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_bbcells_pole_dble.F90 2661 2010-10-25 08:39:13Z coquart $
!  $Author: coquart $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_bbcells_pole_dble.F90 2661 2010-10-25 08:39:13Z coquart $'
!
!----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9990, trim(ch_id)

      call psmile_flushstd
#endif /* VERBOSE */

!
!  Initialization
!
!  MoHa: do I have to take z into account?
      ierror = 0
      method_pole_array_size = 0
      method_pole_array_allocated_size = size(pole_array)
      allocate(method_pole_array_x(method_pole_array_allocated_size))
      allocate(method_pole_array_y(method_pole_array_allocated_size))
      method_pole_array_x = prism_undefined
      method_pole_array_y = prism_undefined
!
!  Detect pole cells, the pole_array contains the corner cells which cover the pole,
!  but we are working with method cells, which are related to the corners cells but not identical
!
#ifdef DEBUG
      PRINT *, TRIM(ch_id), 'number of points at pole found in set_corners :', &
               SIZE(pole_array)
      call psmile_flushstd
#endif
      do n = 1, size(pole_array)

         ijk = psmile_transform_index_1d_to_3d (pole_array(n), corner_shape)

         do i = max (ijk(1)-1, coords_shape(1,1)), min (ijk(1)+1, coords_shape(2,1)-1)

            do j = max (ijk(2)-1, coords_shape(1,2)), min (ijk(2)+1, coords_shape(2,2)-1)
               cell(1)%x = coords_x(i  , j  )
               cell(2)%x = coords_x(i+1, j  )
               cell(3)%x = coords_x(i+1, j+1)
               cell(4)%x = coords_x(i  , j+1)

               cell(1)%y = coords_y(i  , j  )
               cell(2)%y = coords_y(i+1, j  )
               cell(3)%y = coords_y(i+1, j+1)
               cell(4)%y = coords_y(i  , j+1)

               call psmile_transform_cell_cyclic ( cell(:)%x, period, ierror )

               if (psmile_cell_covers_pole(cell(:)%x, cell(:)%y, 4)) then

                  !if the cell is already in the pole array
                  if (.not. any((method_pole_array_x == i) .and. (method_pole_array_y == j))) then

                     if (method_pole_array_size == method_pole_array_allocated_size) then

                        method_pole_array_allocated_size = method_pole_array_allocated_size + 8
                        method_pole_array_x => psmile_realloc(method_pole_array_x, &
                                                            method_pole_array_allocated_size)
                        method_pole_array_y => psmile_realloc(method_pole_array_y, &
                                                            method_pole_array_allocated_size)
                     endif

#if defined(VERBOSE) && defined(DEBUG)
                     print *, trim(ch_id), ": found pole method cell at", i, j
                     print "(a, ': x', /, 2f8.2, /, 2f8.2, /, 2f8.2, /, 2f8.2)", trim(ch_id)
                     PRINT*, 'cell x :',cell(1)%x,cell(2)%x,cell(3)%x,cell(4)%x
                     PRINT*, 'cell y :',cell(1)%y,cell(2)%y,cell(3)%y,cell(4)%y
#endif

                     method_pole_array_size = method_pole_array_size + 1
                     method_pole_array_x(method_pole_array_size) = i
                     method_pole_array_y(method_pole_array_size) = j
                  endif
               endif
            enddo
         enddo
      enddo
!
!  Main part
!
      do n = 1, method_pole_array_size

         i = method_pole_array_x(n)
         j = method_pole_array_y(n)

         ! extending longitude bouding box
         chmax_x(i, j) = 180.0d0
         chmin_x(i, j) = -180.0d0

         ! extent chmin/max to -90 or 90 respectively and adjust midpoint
         if (chmax_y(i, j) > 0) then
             chmax_y(i, j) = 90.0d0

#if defined(VERBOSE) && defined(DEBUG)
         print *, trim(ch_id), ': extending chmax_x to 90'
#endif /* VERBOSE */

         else
            chmin_y(i, j) = -90.0d0

#if defined(VERBOSE) && defined(DEBUG)
         print *, trim(ch_id), ': extending chmin_y to -90'
#endif /* VERBOSE */
         endif

         cell(1)%x = coords_x(i  , j  )
         cell(2)%x = coords_x(i+1, j  )
         cell(3)%x = coords_x(i+1, j+1)
         cell(4)%x = coords_x(i  , j+1)

         cell(1)%y = coords_y(i  , j  )
         cell(2)%y = coords_y(i+1, j  )
         cell(3)%y = coords_y(i+1, j+1)
         cell(4)%y = coords_y(i  , j+1)

         call compute_midpoint(cell, &
            midp_x(i, j), midp_y(i, j))

      enddo
!
!  cleaning up
!
   deallocate (method_pole_array_x, method_pole_array_y)
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

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

      contains

         subroutine compute_midpoint(cell, midp_x, midp_y)
            type(point_dble), intent(in)  :: cell(4)
            double precision, intent(out) :: midp_x, midp_y
            type(point_dble)              :: temp_cell(4), temp_midp

#ifdef VERBOSE
            print 9970, trim(ch_id), cell(1)%x,cell(1)%y
            print 9971, trim(ch_id), cell(2)%x,cell(2)%y
            print 9972, trim(ch_id), cell(3)%x,cell(3)%y
            print 9973, trim(ch_id), cell(4)%x,cell(4)%y
            call psmile_flushstd
#endif /* VERBOSE */

            call psmile_transrot_dble(cell(1)%x, cell(1)%y, temp_cell(1)%x, temp_cell(1)%y)
            call psmile_transrot_dble(cell(2)%x, cell(2)%y, temp_cell(2)%x, temp_cell(2)%y)
            call psmile_transrot_dble(cell(3)%x, cell(3)%y, temp_cell(3)%x, temp_cell(3)%y)
            call psmile_transrot_dble(cell(4)%x, cell(4)%y, temp_cell(4)%x, temp_cell(4)%y)

            temp_midp%x = sum(temp_cell(:)%x) / 4.0d0
            temp_midp%y = sum(temp_cell(:)%y) / 4.0d0

            call psmile_transrot_back_dble(midp_x, midp_y, &
                                           temp_midp%x, temp_midp%y)

#ifdef VERBOSE
            print 9960, trim(ch_id), midp_x, midp_y
            call psmile_flushstd
#endif /* VERBOSE */
9970 format (1x, a, ': psmile_bbcells_pole_dble: computing new midpoint for cell(1) =', 2f8.2)
9971 format (1x, a, ': psmile_bbcells_pole_dble: computing new midpoint for cell(2) =', 2f8.2)
9972 format (1x, a, ': psmile_bbcells_pole_dble: computing new midpoint for cell(3) =', 2f8.2)
9973 format (1x, a, ': psmile_bbcells_pole_dble: computing new midpoint for cell(4) =', 2f8.2)
9960 format (1x, a, ': psmile_bbcells_pole_dble: new midpoint =', 2f8.2)

         end subroutine compute_midpoint

         function psmile_cell_covers_pole (corners_x, corners_y, num_corners)
            double precision, intent(in) :: corners_x(4), corners_y(4)
            integer, intent(in) :: num_corners
            logical :: psmile_cell_covers_pole
            double precision :: angle
            integer :: sense, i

            ! tests if any corner is directly on the pole
            if (any(abs(corners_y) == 90)) then
               psmile_cell_covers_pole = .true.
               return
            endif

            sense = 0
            ! for all edges
            do i = 1, num_corners - 1
               angle = corners_x(i+1) - corners_x(i)
               sense = sense + sign(1.0d0,mod(angle-180.0d0,360.0d0)+180.0d0)
            enddo

            angle = corners_x(1) - corners_x(num_corners)
            sense = sense + sign(1.0d0,mod(angle-180.0d0,360.0d0)+180.0d0)

            psmile_cell_covers_pole = abs(sense) == num_corners
         end function psmile_cell_covers_pole

      end subroutine psmile_bbcells_pole_dble