psmile_remove_intersect.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_Remove_Intersect
!
! !INTERFACE:

      subroutine psmile_remove_intersect (inter, idl, idg, npart,  &
                      local_extent_info, global_extent_info, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_Remove_intersect

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)              :: 
         local_extent_info (nd_extent_infos, *)
!
!     Local extents info with
!         local_extent_info (3, :) = Transformation code
!        
!
      Integer, Intent (In)              :: 
         global_extent_info (nd_extent_infos, *)
!
!     Global extent info with
!         global_extent_info (3, :) = Transformation code
!
!
! !INPUT/OUTPUT PARAMETERS:
!
      Integer, Intent (InOut)           :: npart

!     Number of intersections
!
      Real (PSMILe_float_kind), Intent (InOut) :: inter (2, ndim_3d, npart)

!     Intersections
!
      Integer, Intent (InOut)           :: idl (npart)

!     First (local) Id of the intersections
!     (to be changed if a intersection is removed)
!
      Integer, Intent (InOut)           :: idg (npart)

!     Second (global) Id of the intersections
!     (to be changed if a intersection is removed)
!
! !OUTPUT PARAMETERS:
!
      Integer, Intent (Out)             :: ierror

!     Returns the error code of PSMILe_Remove_intersect;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !LOCAL PARAMETERS
!
      Double Precision                  :: gap = 0.0d0
!
! !LOCAL VARIABLES
!
      Integer                           :: i, k, l, lbeg
      Integer                           :: nfit
      Logical                           :: fit (ndim_3d)
!
      Integer                           :: trans_diff (npart)
      Real (PSMILe_float_kind)          :: common (2, ndim_3d, npart)
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_Remove_intersect" removes the common parts
! from the intersections transferred.
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 03.07.06    H. Ritzdorf  created
!
!EOP
!----------------------------------------------------------------------
!
! $Id: psmile_remove_intersect.F90 2325 2010-04-21 15:00:07Z valcke $
! $Author: valcke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_remove_intersect.F90 2325 2010-04-21 15:00:07Z valcke $'
!
!----------------------------------------------------------------------
!
!  Initialization
!
#ifdef VERBOSE
      print 9990, trim(ch_id), npart

      call psmile_flushstd

#endif /* VERBOSE */
!
      ierror = 0
!
         do l = 1, npart
         trans_diff (l) = local_extent_info(3, idl(l)) - &
                         global_extent_info(3, idg(l))
         end do
!
!===> Look for common intersections
!
      k = 1
          do while (k <= npart - 1)
!
!===> ... Get common part of inter (:, :, k) with other inter's
!
             do l = k+1, npart
             common (1, :, l) = max (inter (1,:,l), inter (1,:,k))
             common (2, :, l) = min (inter (2,:,l), inter (2,:,k))
             end do ! l
!
!     ... Ignore parts with different transformations
!         ? Wird dann noch etwas removet; Dann kann man sich die
!           Routine sparen !
!
             do l = k+1, npart
             if (trans_diff (l) /= trans_diff (k)) &
                common (1, 1, l) = common (2, 1, l) + 1.0
             end do ! l
!
          lbeg = k + 1
!
             do while (lbeg <= npart )
                do l = lbeg, npart
                if (minval (common (2, :, l) - common (1, :, l)) >= 0.0d0) exit
                end do
!
             if (l > npart) exit
             lbeg = l + 1
!
             if (maxval (common (2, :, l) - common (1, :, l)) == 0.0d0) cycle
!
!===> ...... common intersection between k and l found
!            (i)   can inter (:,:,l) be removed ?
!            (ii)  can inter (:,:,k) be removed ?
!            (iii) can inter (:,:,l) be shrinked ?
!            (iv)  ???can inter (:,:,k) be shrinked ?
!
                do i = 1, ndim_3d
                fit (i) = maxval (abs (common (1:2,i,l) - inter (1:2,i,l))) &
                                 == 0.0d0
                end do
             nfit = count (fit)
!
!===> ...... (i) nfit = 3: inter (:, :, l) can be removed
!
             if (nfit == ndim_3d) then
                if (l /= npart) then
                   inter  (:, :, l) = inter  (:, :, npart)
                   common (:, :, l) = common (:, :, npart)
                   idl (l) = idl (npart)
                   idg (l) = idg (npart)
                   trans_diff (l) = trans_diff (npart)
                endif

                npart = npart - 1
                lbeg = l

                cycle ! next l
             endif 
!
!===> ...... (ii) is inter (:, :, k) contained in inter (:, :, l)
!
             if (maxval (abs (common (:,:,l) - inter (:,:,k))) == 0.0d0) then
                inter  (:, :, k) = inter  (:, :, npart)
                trans_diff (k) = trans_diff (npart)

                idl (k) = idl (npart)
                idg (k) = idg (npart)

                npart = npart - 1
                k = k - 1

                exit ! next k
             endif
!
!===> ...... (iii) Can inter (:, :, l) be shrinked ?
!                  changed = .true. : inter (:, :, l) was shrinked
!
!            changed = .false.

             if (nfit >= ndim_2d) then
                if (fit(1) .and. fit (2)) then

                   if      (common (1, 3, l) == inter (1, 3, l)) then
                      inter (1, 3, l) = common (2, 3, l) + gap
!                     changed = .true.

                   else if (common (2, 3, l) == inter (2, 3, l)) then

                      inter (2, 3, l) = common (1, 3, l) - gap
!                     changed = .true.
                   endif

                else if (fit (2) .and. fit (3)) then

                   if      (common (1, 1, l) == inter (1, 1, l)) then
                      inter (1, 1, l) = common (2, 1, l) + gap
!                     changed = .true.

                   else if (common (2, 1, l) == inter (2, 1, l)) then

                      inter (2, 1, l) = common (1, 1, l) - gap
!                     changed = .true.
                   endif

                else ! if (fit (1) .and. fit (3)) then

                   if      (common (1, 2, l) == inter (1, 2, l)) then
                      inter (1, 2, l) = common (2, 2, l) + gap
!                     changed = .true.

                   else if (common (2, 1, l) == inter (2, 1, l)) then

                      inter (2, 2, l) = common (1, 2, l) - gap
!                     changed = .true.
                   endif
                endif
             endif

#ifdef WIRKLICH
!
! in this case, the changed = .true. and .false. statements have to 
! be enabled.
!
!
!===> ... Is shrinking in 3rd Dimension possible
!         Currently, there ist any PRISM grid type where this is necessary !
!
             if (! changed) then

                do i = 1, ndim_3d
                fitk (i) = maxval (abs (common (1:2,i,l) - inter (1:2,i,k))) &
                           == 0
                end do
             nfitk = count (fitk)
!
             if (nfitk >= ndim_2d) then
                if (fitk(1) .and. fitk (2)) then

                   if      (common (1, 3, l) == inter (1, 3, k)) then
                      inter (1, 3, k) = common (2, 3, l) + gap
                      k = k - 1
                      exit ! k

                   else if (common (2, 3, l) == inter (2, 3, k)) then

                      inter (2, 3, k) = common (1, 3, l) - gap
                      k = k - 1
                      exit ! k

                   endif

                else if (fitk (2) .and. fitk (3)) then

                   if      (common (1, 1, l) == inter (1, 1, k)) then
                      inter (1, 1, k) = common (2, 1, l) + gap
                      k = k - 1
                      exit ! k

                   else if (common (2, 1, l) == inter (2, 1, k)) then

                      inter (2, 1, k) = common (1, 1, l) - gap
                      k = k - 1
                      exit ! k
                   endif

                else ! if (fitk (1) .and. fitk (3)) then

                   if      (common (1, 2, l) == inter (1, 2, k)) then
                      inter (1, 2, k) = common (2, 2, l) + gap
                      k = k - 1
                      exit ! k

                   else if (common (2, 1, l) == inter (2, 1, k)) then

                      inter (2, 2, k) = common (1, 2, l) - gap
                      k = k - 1
                      exit ! k
                   endif
                endif
             endif
             endif
#endif
!
             end do ! while (lbeg <= npart)
!
         k = k + 1
         end do ! while ( k <= npart )
!
!===> All done
!
#ifdef VERBOSE
      print 9980, trim(ch_id), ierror, npart
      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!

#ifdef VERBOSE

9990 format (1x, a, ': psmile_remove_intersect: npart =', i4)
9980 format (1x, a, ': psmile_remove_intersect: eof ierror =', i3, &
                    '; npart =', i4)

#endif /* VERBOSE */ 

      end subroutine PSMILe_Remove_intersect