psmile_control_cell_gauss2_dble.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_Control_cell_gauss2_dble
!
! !INTERFACE:

      subroutine psmile_control_cell_gauss2_dble (grid_id,           &
                      ic, nc, list, j, k,                            &
                      coords1,  coords2,  shape,                     &
                      corners1, corners2, corner_shape, nbr_corners, &
                      tol, fnd, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Control_cell_gauss2_dble
#ifdef DEBUG_TRACE
      use psmile_debug_trace
#endif

      Implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: grid_id
!
!     Grid Id (for info degraded cells)
!
      Integer, Intent (In)            :: nc
!
!     Number of cells (and points) to be controlled
!
      Integer, Intent (In)            :: ic (nc)
!
!     Indices of the (corner) cells to be controlled
!     Note: These indices are GaussReduced indices
!
      Integer, Intent (In)            :: list (nc)
!
!     I-Indices of the points in to be controlled
!
      Integer, Intent (In)            :: j
!
!     J-Index of the points in to be controlled
!
      Integer, Intent (In)            :: k
!
!     K-Index of the points in to be controlled
!     Note: list(:), j, k are indices of points to be searched !
!
      Integer, Intent (In)            :: shape (2, ndim_3d)

!     Dimension of coordinate arrays "coord1", "coord2"
!     which contain the coordinates to be searched.

      Double Precision, Intent (In)   ::   
         coords1 (shape(1,1):shape(2,1), shape(1,2):shape(2,2), 
                  shape(1,3):shape(2,3))

      Double Precision, Intent (In)   ::   
         coords2 (shape(1,1):shape(2,1), shape(1,2):shape(2,2), 
                  shape(1,3):shape(2,3))

!     Coordinates to be searched

      Integer, Intent (In)            :: corner_shape (2, ndim_2d)

!     Dimension of corner coordinate arrays

      Integer, Intent (In)            :: nbr_corners

!     Number of corners
!     Note: This routine assumes nbr_corners == 2.

      Double Precision,   Intent (In) ::                 
         corners1 ( corner_shape(1,1):corner_shape(2,1), nbr_corners)

!     Coordinates of corner points of control volume (Longitude)

      Double Precision,   Intent (In) ::                 
         corners2 ( corner_shape(1,1):corner_shape(2,1), nbr_corners)

!     Coordinates of corner points of control volume (Latitude)

      Double Precision, Intent (In)   :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0
!
! !OUTPUT PARAMETERS:
!
      Logical, Intent (Out)           :: fnd (nc)
!
!     Was the i-th point of "list" found within cell "ic"
!
      Integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILe_Control_cell_gauss2_dble;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NTET = Number of triangles used to determine whether
!         point is within the grid cell
!  EPS  = Tolerance for determinate of deformed grid cells
!
      Integer, Parameter              :: ntet = 3
!
      Double Precision, Parameter     :: eps = 1.0d-9
      Double Precision, Parameter     :: rmxeps = 1.0d30
      Double Precision, Parameter     :: dzero = 0.0
!
! !LOCAL VARIABLES
!
!     ... for tolerances
!
      Double Precision                :: tol1
!
!     ... for locations controlled
!
      Integer                         :: i, jj, m, n
!
!     ... Logical reference point and edges (constant)
! 
      Integer, Save                   :: ijkedg (2, ntet)
!
!     ... Reference point and edges and distances
! 
      Double Precision                :: xyzpnt (nc, ntet, ndim_2d), 
                                         xyzref (nc, ntet+1, ndim_2d)
      Double Precision                :: xyzedg (nc, ntet, ndim_2d, 2)
      Double Precision                :: d1, d2
      Double Precision                :: det (nc, ntet)
      Double Precision                :: dnorm
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Control_cell_gauss2_dble" controls whether the point
! "(list(i),j,k)" is contained in the corner cell "ic(i)".
! The routine assumes that the number of corners for the 2-dimensional
! cell is 4 (2 in longitude and 2 in latitude direction).
!
! Subroutine "PSMILe_Control_cell_gauss2_dble" is designed for grids of
! type "PRISM_Gaussreduced_regvrt".
!
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 27.05.08    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_control_cell_gauss2_dble.F90 2927 2011-01-28 14:04:12Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_control_cell_gauss2_dble.F90 2927 2011-01-28 14:04:12Z hanke $'
!
!----------------------------------------------------------------------
!
!  Relative indices of triangles
!
!     The cell is split into 3 triangles
!
!                   4-------3
!                   !       !
!                   !       !
!            Corner 1 ----- 2
!
      data ((ijkedg (jj, n), jj = 1, 2), n = 1, ntet) &
           / 2, 4,  3, 1,  4, 2/
!   Corner    1      2      3
!   Edges of corner 1 are connections to Corner 2 and Corner 4
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), nc

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (tol < 0.0) then
         call psmile_assert (__FILE__, __LINE__, &
                             "tol must be >= 0.0")
      endif
!
!cdir vector
         do i = 1, nc
         if (ic(i) < corner_shape(1,1) .or. ic(i) > corner_shape(2,1) ) then
            exit
         endif
         end do
!
      if (i <= nc) then
         print *, 'i, ic', i, ic(i)
         print *, 'corner_shape', corner_shape
         call psmile_assert (__FILE__, __LINE__, &
                             'wrong index')
      endif
#endif
!
      ierror = 0
!
!===> Compute reference points, corners and edges of transformed
!     triangles.
!
      xyzref (:, 1, 1) = corners1 (ic(:), 1)
      xyzref (:, 1, 2) = corners2 (ic(:), 1)
!
      xyzref (:, 2, 1) = corners1 (ic(:), 2)
      xyzref (:, 2, 2) = xyzref (:, 1, 2)
!
      xyzref (:, 3, 1) = xyzref (:, 2, 1)
      xyzref (:, 3, 2) = corners2 (ic(:), 2)
!
      xyzref (:, 4, 1) = xyzref (:, 1, 1)
      xyzref (:, 4, 2) = xyzref (:, 3, 2)
!
#ifdef DEBUG_TRACE
      do i = 1, nc
        if ( list(i) == ictl_ind(1) .and. j == ictl_ind(2) ) then
           print *, 'Search coords ', &
               coords1(list(i),j,shape(1,3)), coords2(list(i),j,shape(1,3))
           print *, 'Corners: '
           print *, '    C1 : ', xyzref (i, 1, 1), xyzref (i, 1, 2)
           print *, '    C2 : ', xyzref (i, 2, 1), xyzref (i, 2, 2)
           print *, '    C3 : ', xyzref (i, 3, 1), xyzref (i, 3, 2)
           print *, '    C4 : ', xyzref (i, 4, 1), xyzref (i, 4, 2)
        endif
      enddo
#endif
        do n = 1, ntet
        xyzpnt (:, n, 1) = coords1 (list(:), j, k) - xyzref (:, n, 1)
        xyzpnt (:, n, 2) = coords2 (list(:), j, k) - xyzref (:, n, 2)
!
           do m = 1, 2
!cdir vector
              do i = 1, nc
              xyzedg (i, n, 1, m) = xyzref (i, ijkedg(m, n), 1) &
                                  - xyzref (i, n, 1)

              xyzedg (i, n, 2, m) = xyzref (i, ijkedg(m, n), 2) &
                                  - xyzref (i, n, 2)
              end do
           end do
        end do ! n
!
!===> Solve bilinear system; compute spat product
!
!   (xyzedg(*,1), xyzedg(*,2)) (x) = xyzpnt
!
      det(:, :) = xyzedg(:, :, 1,1)*xyzedg(:, :, 2,2) &
                - xyzedg(:, :, 1,2)*xyzedg(:, :, 2,1)
!
#define ALLOW_DEGRADED_CELL
#ifdef ALLOW_DEGRADED_CELL
!
! ??? ist das koscher ?? wird d1(n), d2(n), d3(n) == 0
!
         do n = 1, ntet
!cdir vector
            do i = 1, nc
            dnorm = abs(xyzedg(i, n,1,1)) + abs(xyzedg(i, n,2,1)) + &
                    abs(xyzedg(i, n,1,2)) + abs(xyzedg(i, n,2,2))

            if (abs(det(i, n)) < eps*dnorm*dnorm .or. dnorm == dzero) then
               det (i, n) = sign (rmxeps, det(i, n))
            else
               det (i, n) = 1.0d0 / det (i, n)
            endif
            end do ! i
         end do ! n

#else 

loop_n:  do n = 1, ntet
!cdir vector
            do i = 1, nc
            dnorm = abs(xyzedg(i, n,1,1)) + abs(xyzedg(i, n,2,1)) + &
                    abs(xyzedg(i, n,1,2)) + abs(xyzedg(i, n,2,2))

            if (abs(det(i, n)) < eps*dnorm*dnorm .or. &
                dnorm(i,n) == dzero) then
               det (i, n) = dnorm * 1e-25
! eigentlich exit loop_n und print_degraded cell
! ist bei vektorisierung verschwunden
! sollte man in part_list suchen
            endif
            end do ! i
         end do ! n
!
#ifdef PRINT_DEGRADED_CELL
      if (n < ntet) then
!
!===> cell is degraded
!
         print 9960, trim(ch_id), grid_id, ic (i), &
               (xyzref (i, n, 1), xyzref (i, n, 2), n = 1, 4)
!
            do n = 1, ntet
            if (abs(det(i, n)) < eps) then
               print 9950, n, det (i, n), dnorm, &
                 xyzedg (i, n,1, 1), xyzedg (i, n,1, 2), &
                 xyzedg (i, n,2, 1), xyzedg (i, n,2, 2)
            endif
            end do
!
         nicht gefunden fnd (i) = .false. versucht, durch
         ganz kleine Determinante zu loesen
      endif
#endif /* PRINT_DEGRADED_CELL */
!
      det (:, :) = 1.0d0 / det (:, :)
!
#endif /* ALLOW_DEGRADED_CELL */
!
!===> Return mask for points "list(i)"
!     which were found in the specific cell (ic(i,*), j,k)
!
      tol1 = tol + 1.0
      fnd (:) = .false.
!
         do n = 1, ntet
!cdir vector
            do i = 1, nc

            d1 = (xyzpnt(i, n,1)*xyzedg(i, n,2,2)- &
                  xyzpnt(i, n,2)*xyzedg(i, n,1,2)) * det (i, n)
            d2 = (xyzpnt(i, n,2)*xyzedg(i, n,1,1)- &
                  xyzpnt(i, n,1)*xyzedg(i, n,2,1)) * det (i, n)
!
            fnd (i) = fnd (i) .or. &
               (min(d1, d2) >= -tol .and. &
                max(d1,dzero) + max(d2,dzero) <= tol1)
            end do ! nc
         end do ! ntet
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), count(fnd(:))

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

9990 format (1x, a, ': psmile_control_cell_gauss2_dble: nc', i5)
9980 format (1x, a, ': psmile_control_cell_gauss2_dble: eof, suc', i5)

#ifdef PRINT_DEGRADED_CELL
9960 format (1x, a, ': psmile_control_cell_gauss2_dble: grid id =', i5, &
            /1x,  ' Grid cell (', i4,
                 ') is degraded! Cell coordinates are:', 1p,      &
            4 (/1x, '(', 1 (e17.9, ','), e17.9, ')'))
9950 format (1x, i1, '-th spawn; det =', 1p, e17.9, ' (', e17.9, &
             ') of', 2 (/1x, 2e17.9, ','))
#endif

      end subroutine psmile_control_cell_gauss2_dble