psmile_mg_cells_1d_real.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_MG_cells_1d_real
!
! !INTERFACE:

      subroutine psmile_mg_cells_1d_real (                &
                           nlev, found, loc, range,       &
                           search_grid_type,              &
                           corners1, search_dim,          &
                           shape, control,                &
                           grid_valid_shape, cyclic,      &
                           chmin, chmax,                  &
                           tol, ierror)
!
! !USES:
!
      use PRISM_constants
!
      use PSMILe, dummy_interface => PSMILe_mg_cells_1d_real

      implicit none
!
! !INPUT PARAMETERS:
!
      Integer, Intent (In)            :: nlev

!     Total number of levels

      Integer, Intent (In)            :: search_grid_type

!     Type of source grid

      Integer, Intent (In)            :: range (2, ndim_3d)

!     Dimension of loc and found

      Integer, Intent (In)            :: shape (2, ndim_3d)

!     Dimension of coordinate arrays

      Integer, Intent (InOut)         :: found (range(1,1):range(2,1), 
                                                range(1,2):range(2,2), 
                                                range(1,3):range(2,3))

!     Finest level number on which a grid cell was found for point i,j,k.
!     Input :
!     Level number < -1: Point was not found and
!                        and last level number was (-found(i,j,k))
!     Level number = -(nlev+1): Never found
!     Output :
!     found(i,j,k) = +1: Point (i,j,k) is located on a point
!
!     found(i,j,k) = -1: Point (i,j,k) is located in a cell
!
!     Otherwise, not contained in the cell grid.

      Integer, Intent (InOut)         :: loc   (range(1,1):range(2,1), 
                                                range(1,2):range(2,2), 
                                                range(1,3):range(2,3))

!     Indices of the grid cell in which the point was found.
!     The indices are relative to "shape".

      Real, Intent (In)               :: corners1 (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)            :: control (2, ndim_3d)

!     Input  value: index range in "coords" to be searched
!     Output value: index range in "coords" for which a cell was found

      Integer,          Intent (In)   :: search_dim

!     Dimension in which to be searched

      Integer,          Intent (In)   :: grid_valid_shape (2)

!     Specifies the valid block shape for the "1d"-dimensional block
!     Dimension of chmin and chmax

      Logical,          Intent (In)   :: cyclic

!     Is this a (locally) cyclic coordinate

      Real, Intent (In)               :: chmin (grid_valid_shape(1): 
                                                grid_valid_shape(2)+2)
!
!     Minimum value of bounding box of coordinate.
! chmin (grid_valid_shape(1)+i) contains the minimal value of bounding box of
!                   cell (x_corners(grid_valid_shape(1)+i)  :
!                         x_corners(grid_valid_shape(1)+i+1))
!
! The first box
! chmin (grid_valid_shape(1)-1) contains the minimal value of
!                                        the bounding box from
!                    the virtual cell (x_corners(grid_valid_shape(1)):
!                                               bounding box of
!                                               corner control volume)
! The last box
! chmin (grid_valid_shape(2)) contains the minimal value of
!                                      the bounding box from
!                    the virtual cell (x_corners(grid_valid_shape(2)):
!                                               bounding box of
!                                               corner control volume)
!
      Real, Intent (In)               :: chmax (grid_valid_shape(1): 
                                                grid_valid_shape(2)+2)
!
!     Maximum value of bounding box of coordinate.
!
      Real, Intent (In)               :: tol
!
!     Absolute tolerance for search of "identical" points
!     TOL >= 0.0
!
! !OUTPUT PARAMETERS:
!
      integer, Intent (Out)           :: ierror

!     Returns the error code of PSMILE_mg_cells_1d_real;
!             ierror = 0 : No error
!             ierror > 0 : Severe error
!
! !DEFINED PARAMETERS:
!
!  NREFD = Number of method coordinates to be controlled.
!  val_direct  = Code for locations which should to be directly transferred
!                to the destination process.
!  val_coupler = Code for locations which should to be transferred
!                to the coupler process.

      Integer, Parameter              :: nrefd = 3
!
      Integer, Parameter              :: val_direct  =  1
      Integer, Parameter              :: val_coupler = -1
!
!     ... for levels
!
      Integer, Parameter              :: lev = 1
!
! !LOCAL VARIABLES
!
!     ... for locations searched
!
      Integer                         :: i, j, k, jj
      Integer                         :: ibegl, iendl
      Integer                         :: idim
      Integer                         :: ic
      Integer                         :: ijk(3)
!
!     ... For locations controlled
!
      Integer                         :: irange (2)
      Integer                         :: ctrl(2, ndim_3d)
!
!     ... for search
!
      Integer                         :: ii
!
#ifdef DEBUG
!     ... for locations searched
!
      Integer                         :: nfnd (0:3)
#endif /* DEBUG */
!
! !DESCRIPTION:
!
! Subroutine "PSMILe_mg_cells_1d_real" searches the corresponding cells
! of the grid for the subgrid corners of the sending process.
!
! TODO: Improve !!
!       Use tol(erance) for overlap tests. Same tol has to be used later in
!       PSMILe_Neigh_cells_*
!
! !REVISION HISTORY:
!
!   Date      Programmer   Description
! ----------  ----------   -----------
! 06.07.21    R. Redler    created
!
!EOP
!----------------------------------------------------------------------
!
!  $Id: psmile_mg_cells_1d_real.F90 2082 2009-10-21 13:31:19Z hanke $
!  $Author: hanke $
!
   Character(len=len_cvs_string), save :: mycvs = 
       '$Id: psmile_mg_cells_1d_real.F90 2082 2009-10-21 13:31:19Z hanke $'
!----------------------------------------------------------------------
!
!  Initialization
!
      ierror = 0
      ctrl = control

#ifdef VERBOSE
      print 9990, trim(ch_id), lev, control

      call psmile_flushstd
#endif /* VERBOSE */
!
#ifdef PRISM_ASSERTION
      if (tol < 0.0) then
         call psmile_assert (__FILE__, __LINE__, &
                             "tol must be >= 0.0")
      endif
#endif
!
! create temporary memory (control (2,1)-control(1,1)+1)*ndtry
! and do n =1,nref before i-loop ?!?
!
! ? instead of controlling all nref, determine location in cell once
!
! ? one-dimensional Loop instead of i,j loop
!
#ifdef DEBUG
      nfnd (:) = 0
#endif /* DEBUG */

#ifdef DEBUGX
      print *, 'begin ---'
      do k = control (1, 3), control (2,3)   
         do j = control(1,2), control (2,2) 
            do i = control(1,1), control (2,1) 
               print *, 'i, j, k, corner, loc, found', i,j,k, corners1 (i,j,k), &
                    loc(i,j,k), found (i,j,k)
            end do
         end do
      end do
#endif
      if ( search_grid_type == PRISM_Gaussreduced_regvrt ) then
         idim = 1
      else
         idim = search_dim
      endif
      ijk = 0
      ijk(idim) = 1
      ctrl (2,idim) = control (2,idim) - 1
!
!===>
!
      do k = ctrl (1, 3), ctrl (2,3)
         do j = ctrl(1,2), ctrl (2,2)
!
!-----------------------------------------------------------------------
!        Look in the fine parts of the coarse grid cell loc(i)
!-----------------------------------------------------------------------
!
            ibegl = ctrl (1, 1)
!
!===> ... Look for the next cell to be controlled
!
            do while (ibegl <= ctrl(2,1))
!
!cdir vector
               do i = ibegl, ctrl(2,1)
                  if (found(i, j, k) /= lev) exit
               end do
!
            iendl = i - 1
!
!===> ... Look in the real cell (ic)
!
            do i = ibegl, iendl
!
               ic = loc (i,j,k)
!
               irange (1) = max (ic-1, grid_valid_shape (1))
               irange (2) = min (ic+1, grid_valid_shape (2))

               do ii = irange (1), irange (2)
                  if ( chmin(ii) <= corners1(i,j,k) .and. &
                       chmax(ii) >= corners1(i,j,k) ) exit
               end do

               if ( chmin(ii) > corners1(i,j,k) ) ii = ii - 1

               if (ii <= irange(2)) then

                  loc (i,j,k) = ii
#ifdef DEBUG
                  nfnd (1) = 0
                  nfnd (2) = nfnd (2) + 1
#endif /* DEBUG */
                  found (i,j,k) = val_coupler

               else
!
!              Cell has no overlap with any other cell in the search grid,
!              mark point as 0
#ifdef DEBUG
                  nfnd (3) = nfnd (3) + 1
!                 nfnd (2) = nfnd (2) - 1
#endif /* DEBUG */
!Todo             found (i,j,k) = 0

               endif

            end do ! i
!
!===> ... Is there an cyclic index
!
            if (cyclic) then
!
!              This is an cyclic index;
!              Move points, which are located in virtual cell 0
!              into cell "grid_valid_shape(2)"
!
!cdir vector
               do i = ibegl, iendl
                  if ( loc (i,j,k) < grid_valid_shape (1) ) then
                       loc (i,j,k) = grid_valid_shape (2)
                  endif
               enddo
!
            endif
!
!===> ... Start index of the next search
!
            ibegl = iendl + 2
            end do ! while
!
         end do ! j
      end do ! k
!
!===> Additional check for empty cells. The initial search with mg_method only searched
!     for the minimum corners. Here we check whether the maximum falls into the valid
!     area and in this case mark the cell as found. 
!
      if ( search_grid_type == PRISM_Gaussreduced_regvrt ) then

         jj = (range(2,1)-range(1,1)+1)/2
         ibegl = range(1,1)
         iendl = range(1,1) + jj - 1

         do k = ctrl (1, 3), ctrl (2,3)
            do j = ctrl(1,2), ctrl (2,2)

               do i = ibegl, iendl

                  if ( abs(found(i   ,j,k)) /= 1 .and. &
                       abs(found(i+jj,j,k)) == 1 ) then
                     loc  (i,j,k) = loc  (i+jj,j,k)
                     found(i,j,k) = found(i+jj,j,k)
                  endif

               enddo

               do i = iendl+1, range(2,1)
                  found(i,j,k) = -(nlev+1)
               enddo

            enddo
         enddo

      else

         do k = ctrl (1,3), ctrl (2,3) 
            do j = ctrl (1,2), ctrl (2,2)
               do i = ctrl (1,1), ctrl (2,1)
                  if ( abs(found(i,j,k)) == 1 ) found(i,j,k) = val_coupler
                  if ( abs(found(i,j,k)) > 1 .and. abs(found(i+ijk(1),j+ijk(2),k+ijk(3))) == 1 ) then
                          found(i,j,k) = found(i+ijk(1),j+ijk(2),k+ijk(3))
                          loc(i,j,K)   = loc  (i+ijk(1),j+ijk(2),k+ijk(3))
                  endif
               enddo
            enddo
         enddo

      endif
!
!===> Check that location found are in the valid range, otherwise mark point as not found.
!     Enforce to send all data through the coupler.
!
      do k = range (1,3), range (2,3)   
         do j = range(1,2), range (2,2) 
            do i = range(1,1), range (2,1) 
               if ( abs(found(i,j,k)) == 1 ) found(i,j,k) = val_coupler
               if ( loc(i,j,k) < grid_valid_shape(1) .or. &
                    loc(i,j,k) > grid_valid_shape(2) )    &
                    found(i,j,k) = -(nlev+1)
            end do
         end do
      end do
!
!===> All done
!

#ifdef DEBUGX
      print *, 'end ---'
        do k = control (1, 3), control (2,3)
            do j = control(1,2), control (2,2)
               do i = control(1,1), control (2,1)
               print *, 'i, j, k, corner, loc, found', i,j,k, corners1 (i,j,k), &
                        loc(i,j,k), found (i,j,k)
               end do
            end do
        end do
#endif

#ifdef DEBUG
      print 9970, trim(ch_id), lev, nfnd
9970  format (1x, a, ': psmile_mg_cells_1d_real: lev =', i3, &
              ': not :', i5, ', dir: ', i6, ', fnd :', i6, i5)
#endif /* DEBUG */

#ifdef VERBOSE
      print 9980, trim(ch_id), lev

      call psmile_flushstd
#endif /* VERBOSE */
!
!  Formats:
!
9990 format (1x, a, ': psmile_mg_cells_1d_real: level', i3, &
                    ', control', 6i6)
9980 format (1x, a, ': psmile_mg_cells_1d_real: eof, level', i3)

      end subroutine psmile_mg_cells_1d_real