psmile_celltest_real.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
! Copyright 2007-2010, NEC Europe Ltd., London, UK.
! All rights reserved. Use is subject to OASIS4 license terms.
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: PSMILe_Celltest_real
!
! !INTERFACE:
!
      subroutine psmile_celltest_real ( grid_id, range, sense, ierror )
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Celltest_real

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)          :: grid_id
      Integer, Intent (In)          :: range(2,ndim_2d)
      Integer, Intent (Out)         :: sense
      Integer, Intent (Out)         :: ierror

!     Returns the error code of PSMILe_Celltest_real
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL VARIABLES
!
      Type (Corner_Block), Pointer  :: cp

      Integer                       :: i, j, index
      Integer                       :: ibeg, iend
      Integer                       :: jbeg, jend

      Integer                       :: n, np1, stride, isize
      Integer                       :: nbr_edges

      Integer                       :: sense_p
      Integer                       :: sense_n

      Real, Parameter               :: real_one  = 1.0
      Real, Parameter               :: real_zero = 0.0
      Real, Allocatable             :: edge(:,:)
      Real                          :: cross_prod
!
      Integer, Parameter            :: nerrp = 3
      Integer                       :: ierrp (nerrp)
!
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Celltest_real" checks whether the
! cell information is consistent with PSMILe internal
! assumptions. First we test the cell for convexity.
! Later we can expand this test for valid non-convex
! elements by evaluating the angles of non-convex elements.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 05.11.19    R. Redler    created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_celltest_real.F90 2325 2010-04-21 15:00:07Z valcke $
!  $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_celltest_real.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
#ifdef VERBOSE
      print 9990, trim(ch_id), range
      call psmile_flushstd
#endif /* VERBOSE */
!
!  Initialization
!
      ierror  = 0
      sense_p = 0
      sense_n = 0

      if ( Grids(grid_id)%grid_type /= PRISM_Irrlonlat_regvrt ) then
         ierrp (1) = Grids(grid_id)%grid_type
         ierror = PRISM_Error_Internal
         call psmile_error ( ierror, 'unsupported grid type', &
              ierrp, 1, __FILE__, __LINE__ )
         return
      endif

      cp => Grids(grid_id)%corner_pointer

      ibeg = range(1,1) - cp%corner_shape(1,1) + 1
      iend = range(2,1) - cp%corner_shape(1,1) + 1
      jbeg = range(1,2) - cp%corner_shape(1,2) + 1
      jend = range(2,2) - cp%corner_shape(1,2) + 1

      nbr_edges = cp%nbr_corners/2
      
      stride = (cp%corner_shape(2,1)-cp%corner_shape(1,1)+1) * &
               (cp%corner_shape(2,2)-cp%corner_shape(1,2)+1)

      isize  =  cp%corner_shape(2,1)-cp%corner_shape(1,1)+1
      !
      !----------------------------------
      ! Test for convexity of a the cells
      !   still needs to be vectorised.
      !----------------------------------
      !

      Allocate ( edge(2,nbr_edges) )

      do j = jbeg, jend
      do i = ibeg, iend

         index = (j-1)*isize+i

         do n = 1, nbr_edges-1
            edge(1,n) = cp%corners_real(1)%vector( n   *stride+index) - &
                        cp%corners_real(1)%vector((n-1)*stride+index)
            edge(2,n) = cp%corners_real(2)%vector( n   *stride+index) - &
                        cp%corners_real(2)%vector((n-1)*stride+index)
         enddo

         edge(1,nbr_edges) = cp%corners_real(1)%vector(index) - &
                             cp%corners_real(1)%vector((nbr_edges-1)*stride+index)
         edge(2,nbr_edges) = cp%corners_real(2)%vector(index) - &
                             cp%corners_real(2)%vector((nbr_edges-1)*stride+index)
         !
         !------------------------------------------------------
         ! get sign of cross products of egde vectors and sum up
         !------------------------------------------------------
         !
         sense    = 0

         do n = 1, nbr_edges
            np1 = mod(n, nbr_edges) + 1
            cross_prod = edge(1,n)*edge(2,np1) - edge(2,n)*edge(1,np1)
            if ( cross_prod /= real_zero ) &
            sense = sense + sign(real_one, cross_prod)
         enddo

         if ( abs(sense) /= nbr_edges ) then

            ierrp(1) = i + cp%corner_shape(1,1) - 1
            ierrp(2) = j + cp%corner_shape(1,2) - 1

            ierror = PRISM_Warn_Cell
            call psmile_warning ( ierror, 'Non-convex or degraded cell', &
                 ierrp, 2, __FILE__, __LINE__ )
            ierror = 0
#ifdef DEBUG
            do n = 0, nbr_edges-1
             print *, ' corner: n, lon, lat', n+1, cp%corners_real(1)%vector(n*stride+index),  &
                                                   cp%corners_real(2)%vector(n*stride+index)
            enddo
            print *, ' calculated sense is ', sense
#endif
         endif

         if ( sign(1,sense) > 0 ) then
            sense_p = sense_p + 1
         else
            sense_n = sense_n + 1
         endif

      enddo
      enddo

      Deallocate ( edge )

      if ( sense_p * sense_n /= 0 ) then
         ierror = PRISM_Warn_Cell
         ierrp(1) = sense_p
         ierrp(2) = sense_n
         call psmile_warning ( ierror, 'Sense of cells must be the same', &
            ierrp, 2, __FILE__, __LINE__ )
         ierror = 0
      endif

      sense = sign(1,sense)
!
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990 format (1x, a, ': PSMILe_Celltest_real: range ', 4i8)
9980 format (1x, a, ': PSMILe_Celltest_real: eof ierror =', i3)
!
      end subroutine PSMILe_Celltest_real