psmile_neigh_extra_points.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_Neigh_extra_points
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_extra_points (search, extra_search, &
00012                mask_array, mask_shape, mask_available, use_how,   &
00013                grid_valid_shape,                                  &
00014                neighbors_3d, nloc, n_corners, len_cpl, ierror)
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019 !
00020       use PSMILe, dummy_interface => PSMILe_Neigh_extra_points
00021 #ifdef DEBUG_TRACE
00022       use psmile_debug_trace
00023 #endif
00024 
00025       Implicit none
00026 !
00027 ! !INPUT PARAMETERS:
00028 !
00029       Type (Enddef_search), Intent (In) :: search
00030 
00031 !     Info's on coordinates to be searched
00032 !
00033       Integer,         Intent (In)    :: mask_shape (2, ndim_3d)
00034 
00035 !     Dimension of (method) mask array
00036 
00037       Logical,         Intent (In)    :: mask_array (mask_shape (1,1): 
00038                                                      mask_shape (2,1), 
00039                                                      mask_shape (1,2): 
00040                                                      mask_shape (2,2), 
00041                                                      mask_shape (1,3): 
00042                                                      mask_shape (2,3))
00043 !     Mask of the method
00044 
00045       Logical,         Intent (In)    :: mask_available
00046 
00047 !     Is mask specified in array "mask_array" ?
00048 !     mask_available = .false. : Mask is not available
00049 !                      .true.  : Mask is     available
00050 
00051       Integer,         Intent (In)    :: use_how(3)
00052 
00053 !     Information about how to apply the mask
00054 !        PSMILE_novalue   : If one neighbor is unvalid reject all.
00055 !        PSMILE_tneighbour: Reject all unvalid neighbors but keep the valid.
00056 !        PSMILE_nneighbour: == PSMILE_undef
00057 
00058       Integer,         Intent (In)    :: grid_valid_shape (2, ndim_3d)
00059 
00060 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00061 !
00062       Integer,         Intent (In)    :: nloc
00063 !
00064 !     Total number of locations to be transferred.
00065 !     Number of locations for which extra points should be searched.
00066 !     Second dimension of neighbors array "neighbors_3d"
00067 !
00068       Integer,         Intent (In)    :: n_corners
00069 !
00070 !     Last dimension of neighbors array "neighbors_3d" and
00071 !     number of interpolation bases.
00072 !
00073       Integer, Intent (In)            :: len_cpl (search%search_data%npart)
00074 
00075 !     Number of points to be sent to the coupler for each partition.
00076 !
00077 ! !INPUT/OUTPUT PARAMETERS:
00078 !
00079       Type (Extra_search_info), Intent (InOut) :: extra_search
00080 !
00081 !     Number of locations where
00082 !       (*) required mask values were not true
00083 !
00084       Integer,         Intent (InOut) :: neighbors_3d (ndim_3d, nloc, n_corners)
00085 
00086 !     Indices of neighbor locations (interpolation bases)
00087 !
00088 ! !OUTPUT PARAMETERS:
00089 
00090       Integer,         Intent (Out)   :: ierror
00091 
00092 !     Returns the error code of PSMILE_Neigh_extra_points;
00093 !             ierror = 0 : No error
00094 !             ierror > 0 : Severe error
00095 !
00096 ! !LOCAL VARIABLES
00097 !
00098 !     ... for locations searched
00099 !
00100       Integer                         :: grid_id
00101       Integer                         :: search_mode
00102       Integer                         :: i, j, n, nbeg, nend
00103 !
00104 !     ... for locations controlled
00105 !
00106       Integer                         :: n_extra, outside
00107 !
00108       Integer, Pointer                :: indices (:)
00109       Logical, Allocatable            :: mask_cell (:)
00110 !
00111 !     ... for partitions
00112 !
00113       Integer                         :: ipart
00114 !
00115 !     ... for error handling
00116 !
00117       Integer, parameter              :: nerrp = 2
00118       Integer                         :: ierrp (nerrp)
00119 !
00120 !
00121 ! !DESCRIPTION:
00122 !
00123 ! Subroutine "PSMILe_Neigh_extra_points" determines the points (indices)
00124 ! which require extra search (nearest neighbour search).
00125 ! The extra points are searched the neighbour indices (interpolation
00126 ! bases) specified in
00127 !
00128 !      neighbour (1:ndim_3d, 1:nloc, n_corners).
00129 !
00130 ! !REVISION HISTORY:
00131 !
00132 !   Date      Programmer   Description
00133 ! ----------  ----------   -----------
00134 ! 15.03.06    H. Ritzdorf  created
00135 !
00136 !EOP
00137 !----------------------------------------------------------------------
00138 !
00139 !  $Id: psmile_neigh_extra_points.F90 3174 2011-05-04 09:41:43Z hanke $
00140 !  $Author: hanke $
00141 !
00142    Character(len=len_cvs_string), save :: mycvs = 
00143        '$Id: psmile_neigh_extra_points.F90 3174 2011-05-04 09:41:43Z hanke $'
00144 !
00145 !----------------------------------------------------------------------
00146 !
00147 !  Initialization
00148 !
00149 #ifdef VERBOSE
00150       print 9990, trim(ch_id)
00151 
00152       call psmile_flushstd
00153 #endif /* VERBOSE */
00154 !
00155 !===> Initialization
00156 !
00157       ierror  = 0
00158       search_mode = use_how(1)
00159       grid_id = use_how(3)
00160 !
00161 !     Note: "outside" must be greater than
00162 !           "grid_valid_shape (2,1) + maximal with of interpolation stencil"
00163 !            in order to enable global search.
00164 !
00165       if (associated(Grids(grid_id)%reduced_gauss_data%remote_index)) then
00166           !
00167           ! in particular for Gauss-reduced grids
00168           !
00169           outside   = grid_valid_shape(2,1) + &
00170                       size(Grids(grid_id)%reduced_gauss_data%remote_index) + 11
00171       else
00172           outside   = grid_valid_shape(2,1) + 11
00173       endif
00174 !
00175 #ifdef PRISM_ASSERTION
00176 !
00177       if (nloc /= Sum(len_cpl)) then
00178          print *, "nloc, Sum(len_cpl)", nloc, Sum (len_cpl)
00179          call psmile_assert (__FILE__, &
00180               __LINE__, "ncpl == SUM(len_cpl) expected")
00181       endif
00182 
00183       if (extra_search%n_extra /= 0) then
00184          call psmile_assert (__FILE__, &
00185               __LINE__, "extra_search%n_extra /= 0")
00186       endif
00187 #endif
00188 !
00189 !-----------------------------------------------------------------------
00190 !     Control Mask values
00191 !
00192 !     The searched point is located in donor cell
00193 !     (srcloc(:,i), srcloc(:,i)+1).
00194 !-----------------------------------------------------------------------
00195 !
00196 ! PSMILE_novalue  : If one neighbor is unvalid reject all.
00197 !                   This will cause the transformer to perform
00198 !                   bi/trilinear interpolation for those points
00199 !                   only for which all neighbors are valid.
00200 !
00201 ! PSMILE_tneighbour: Reject all unvalid neighbors but keep the valid.
00202 !                    This will cause the transformer to perform
00203 !                    bi/trilinear interpolation for those points
00204 !                    only for which all neighbors are valid.
00205 !                    Those for which less than "n_corners" valid neighbors
00206 !                    have been found will be calculated using weighted
00207 !                    distances.
00208 !
00209 ! PSMILE_nneighbour: == PSMILE_undef [default]
00210 !                    Bi/trilinear interpolation for those points for which
00211 !                    all "n_corners" points are valid.
00212 !                    For those for which less than "n_corners" valid
00213 !                    neighbors were found those neighbors are kept.
00214 !                    For the "n_extra" target points for which 0 valid
00215 !                    neighbors were found we invoke the search for
00216 !                    one single nearest neighbor.
00217 !
00218 ! NOTE: If a interpolation base is outside of the physical boundary,
00219 !       the base points is handled such as the mask would be false.
00220 !
00221 !-----------------------------------------------------------------------
00222 !
00223 !     set search mode PSMILE_nneighbour to default
00224 !
00225 #ifdef DEBUG
00226 8990  format (1x, a, ': psmile_neigh_extra_points: search_mode', i8, &
00227               '; mask_available ', l1)
00228       print 8990, trim(ch_id), search_mode, mask_available
00229 #endif
00230 !
00231       if ( search_mode == PSMILE_nneighbour ) search_mode = PSMILe_Undef
00232 
00233 #ifdef DEBUG_TRACE
00234          if ( ictl > 0 .and. ictl <= nloc ) then
00235             print *, ' ictl neigh 1', neighbors_3d(1,ictl,:)
00236             print *, ' ictl neigh 2', neighbors_3d(2,ictl,:)
00237             print *, ' ictl neigh 3', neighbors_3d(3,ictl,:)
00238          endif
00239 #endif
00240  if (mask_available) then
00241 
00242 !        Ignore corners of (method) cells flagged as false by moving them
00243 !        out of the valid range (see psmile_compact_neighbors_3d)
00244 
00245          do n = 1, n_corners
00246 !cdir vector
00247             do i = 1, nloc
00248 
00249               if (neighbors_3d(1, i, n) /= extra_search%global_marker) then
00250 !
00251 !             search_mode == vneighbour, not implemented !
00252 !
00253 !             If the point was not found in regular domain,
00254 !             mark it by "outside".
00255 !
00256                   if ( neighbors_3d(1, i, n) < grid_valid_shape(1,1) .or. &
00257                        neighbors_3d(1, i, n) > grid_valid_shape(2,1) .or. &
00258                        neighbors_3d(2, i, n) < grid_valid_shape(1,2) .or. &
00259                        neighbors_3d(2, i, n) > grid_valid_shape(2,2) .or. &
00260                        neighbors_3d(3, i, n) < grid_valid_shape(1,3) .or. &
00261                        neighbors_3d(3, i, n) > grid_valid_shape(2,3) ) then
00262 !
00263                        neighbors_3d(1, i, n) = outside
00264 !
00265 !             This corresponds to "search_mode == tneighbour".
00266 !             Only the actual neighour point is moved outside
00267 !             since the mask is .false..
00268 !
00269                   else if ( neighbors_3d(1, i, n) >= mask_shape(1,1) .and. &
00270                             neighbors_3d(1, i, n) <= mask_shape(2,1) .and. &
00271                             neighbors_3d(2, i, n) >= mask_shape(1,2) .and. &
00272                             neighbors_3d(2, i, n) <= mask_shape(2,2) .and. &
00273                             neighbors_3d(3, i, n) >= mask_shape(1,3) .and. &
00274                             neighbors_3d(3, i, n) <= mask_shape(2,3) ) then
00275 #ifdef DEBUG_TRACE
00276                      if ( i == ictl ) then
00277                           print *, ' ictl neigh & mask ', i, n, neighbors_3d(:, i, n), &
00278                                                     mask_array (neighbors_3d(1, i, n), &
00279                                                                 neighbors_3d(2, i, n), &
00280                                                                 neighbors_3d(3, i, n))
00281                      endif
00282 #endif
00283                      if ( .not. mask_array (neighbors_3d(1, i, n),   &
00284                                             neighbors_3d(2, i, n),   &
00285                                             neighbors_3d(3, i, n)) ) &
00286                        neighbors_3d(1, i, n) = outside
00287                   endif
00288 
00289                endif
00290             enddo ! i
00291          enddo ! n
00292 !
00293 !    search_mode == tneighbour is already treated above. Now we check whether
00294 !    additional neighbours have to be marked as not valid (outside) for
00295 !        search_mode == novalue
00296 !        search_mode == nneighbour | undef [default]  
00297 !
00298          if ( search_mode == PSMILE_novalue .or. &
00299               search_mode == PSMILe_undef ) then
00300 !
00301 !===> ... Allocate and initialize work vector "mask_cell" with
00302 !         mask_cell (i) = .true. : Mask of first point is .false.
00303 !                                  ==> First index is out of range
00304 !
00305             Allocate (mask_cell (nloc), STAT = ierror)
00306             if (ierror > 0) then
00307                ierrp (1) = ierror
00308                ierrp (2) = nloc
00309 
00310                ierror = PRISM_Error_Alloc
00311 
00312                call psmile_error ( ierror, 'mask_cell', &
00313                     ierrp, 2, __FILE__, __LINE__ )
00314                return
00315             endif
00316 !
00317 !cdir vector
00318             do i = 1, nloc
00319                mask_cell (i) = neighbors_3d(1, i, 1) == outside
00320             end do
00321 !
00322             if ( search_mode == PSMILE_novalue ) then
00323 !
00324 !           Ignore cells with one or more neighbors flagged as false 
00325 !           by moving all neighbors out of the valid range
00326 !           (see psmile_compact_neighbors_3d)
00327 !
00328 !           mask_cell (i) = .true. : A mask entry of (method) cell is .false.
00329 !
00330                do n = 2, n_corners
00331 !cdir vector
00332                   do i = 1, nloc
00333                      mask_cell (i) = mask_cell (i) .or. &
00334                           neighbors_3d(1, i, n) == outside
00335                   end do
00336                end do
00337 
00338                do n = 1, n_corners
00339                   !cdir vector
00340                   do i = 1, nloc
00341                      if ( mask_cell(i) ) &
00342                         neighbors_3d(1, i, n) = outside
00343                   enddo
00344                enddo
00345 !
00346             else ! if ( search_mode == PSMILe_undef ) then
00347 !
00348 !               mask_cell (i) = .true. : All mask entries of (method) cell
00349 !                                        are .false.
00350 !
00351                do n = 2, n_corners
00352 !cdir vector
00353                   do i = 1, nloc
00354                      mask_cell (i) = mask_cell (i) .and. &
00355                           neighbors_3d(1, i, n) == outside
00356                   end do
00357                end do
00358 !
00359 #ifdef TEST_EXTRA_SEARCH
00360                print *, '### psmile_neigh_extra_points: Extra-Search for testing enabled'
00361                mask_cell(:) = .true.
00362 #endif
00363 !
00364 !===> Store cells/points which require extra search
00365 !     Note: The index in the entire vector "neighbors_3d" is stored
00366 !           and not the relative index in the partition "ipart"
00367 !
00368 !===> ... Get number of points to be searched extra
00369 !
00370                n_extra = count (mask_cell (:))
00371 !
00372                if (n_extra > 0) then
00373 !
00374 !===> ...... Allocate memory required to store the indices and
00375 !            code for the points required.
00376 !            Note: The indices are stored as index in entire (!) 
00377 !                  array "neighbors"
00378 !
00379                   Allocate (indices (n_extra), stat = ierror)
00380                   if (ierror > 0) then
00381                      ierrp (1) = ierror
00382                      ierrp (2) = n_extra
00383 
00384                      ierror = PRISM_Error_Alloc
00385 
00386                      call psmile_error ( ierror, 'extra_search%indices', &
00387                           ierrp, 2, __FILE__, __LINE__ )
00388                      return
00389 endif
00390 !
00391 !===> ...... Generate indices
00392 !
00393 #ifdef USE_PACK
00394                   indices (1: n_extra) = &
00395                        PACK ((/ (i, i=1,nloc) /), mask_cell)
00396 #else
00397                   j = 0
00398 !cdir vector
00399                   do i = 1, nloc
00400                      if (mask_cell (i)) then
00401                         j = j + 1
00402                         indices (j) = i
00403                      endif
00404                   end do
00405 #endif /* USE_PACK */
00406 !
00407 !===> ...... Generate length per partition
00408 !
00409                   nend = 0
00410                   do ipart = 1, search%search_data%npart
00411                      if (len_cpl (ipart) > 0) then
00412                         nbeg = nend + 1
00413                         nend = nend + len_cpl (ipart)
00414 !
00415                         extra_search%len_extra (ipart) = &
00416                              count (mask_cell(nbeg:nend))
00417 !
00418                      endif ! len_cpl > 0
00419                   end do ! ipart
00420 !
00421 #ifdef DEBUGX
00422                   do i = 1, nloc
00423                      if (mask_cell (i)) then
00424                         print *, 'i', i, neighbors_3d(:, i, 1)
00425                         do n = 1, n_corners
00426                            print *, 'ng', neighbors_3d(:, i, n), &
00427                                 mask_array (neighbors_3d(1, i, n), &
00428                                 neighbors_3d(2, i, n), &
00429                                 neighbors_3d(3, i, n))
00430                         end do
00431                      endif
00432                   end do
00433 #endif
00434 !
00435 !===> ...... Save extra search data in structure "extra_search"
00436 !
00437                   extra_search%indices => indices
00438 !
00439                   extra_search%n_extra = n_extra
00440 
00441                endif ! (n_extra > 0)
00442             endif ! if (earch_mode == PSMILE_novalue)
00443 !
00444             Deallocate (mask_cell)
00445             !
00446          endif ! if (search_mode == PSMILE_novalue .or. ...
00447 
00448       endif ! if (mask_available)
00449 !
00450 !===> All done
00451 !
00452 #ifdef VERBOSE
00453       print 9980, trim(ch_id), extra_search%n_extra
00454 
00455 !     print *, extra_search%len_extra (:)
00456       call psmile_flushstd
00457 #endif /* VERBOSE */
00458 !
00459 !  Formats:
00460 !
00461 9990  format (1x, a, ': psmile_neigh_extra_points:')
00462 9980  format (1x, a, ': psmile_neigh_extra_points: eof, n_extra', i6)
00463 
00464     end subroutine psmile_neigh_extra_points