psmile_neigh_cells_3d_real.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 #if 0
00008 #define GRIDPOINT 4417
00009 #define META_RANK 1
00010 #define META_POST
00011 #endif
00012 !
00013 ! !ROUTINE: PSMILe_Neigh_cells_3d_real
00014 !
00015 ! !INTERFACE:
00016       subroutine psmile_neigh_cells_3d_real ( use_how,      &
00017                grid_shape, interpolation_mode, cyclic,      &
00018                corner_shape_3d, nbr_corners,                &
00019                corner_x, corner_y, corner_z,                &
00020                search,  control, tgt_cell, tgt_corners,     &
00021                npoints, srclocs, msklocs, ncpl,             &
00022                num_neigh, nbr_cells, ierror )
00023 !
00024 ! !USES:
00025 !
00026       use PRISM_constants
00027 !
00028       use PSMILe, dummy_interface => PSMILe_Neigh_cells_3d_real
00029 
00030       implicit none
00031 !
00032 ! !INPUT PARAMETERS:
00033 !
00034       Integer, Intent (In)            :: use_how(3)
00035 !
00036 !     container of information  under consideration
00037 !
00038       Integer, Intent (In)            :: grid_shape (2, ndim_3d)
00039 !
00040 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00041 !
00042       Integer, Intent (In)            :: interpolation_mode
00043 !
00044 !     Specifies the interpolation mode with
00045 !        interpolation_mode = 3 : 3D conservative remapping
00046 !                                 (not supported)
00047 !        interpolation_mode = 2 : 2D conservative remapping
00048 !        interpolation_mode = 1 : not supported
00049 !                                 (determines 2 corners)
00050 !
00051       Logical, Intent (In)            :: cyclic (ndim_3d)
00052 !    
00053 !     Does the block have cyclic coordinates in the corresponding direction ?
00054 !
00055       Integer, Intent (In)            :: corner_shape_3d (2,ndim_3d,ndim_3d)
00056 !
00057 !     Shape of local corners -> corner_x, corner_y
00058 !
00059       Integer, Intent (In)            :: nbr_corners(ndim_3d)
00060 !
00061 !     Number of corners of local corners corner_x, corner_y, corner_z
00062 !
00063       Type (Enddef_search), Intent (InOut) :: search
00064 !
00065 !     Info's on coordinates to be searched
00066 !
00067       Integer, Intent (In)            :: control (2, ndim_3d, search%npart)
00068 !
00069 !     Control range
00070 !
00071       Type (real_vector), Intent (InOut) :: tgt_cell(ndim_3d)
00072 !
00073 !     Cells for which we search
00074 !
00075       Integer, Intent (In)            :: tgt_corners(ndim_3d)
00076 !
00077 !     Number of corner in the target cell
00078 !
00079       Integer, Intent (In)            :: npoints
00080 !
00081 !     Total number of locations to be transferred
00082 !
00083       Integer, Intent (In)            :: ncpl
00084 !
00085 !     Total number of locations
00086 !
00087       Integer, Intent (In)            :: srclocs (ndim_3d, npoints)
00088 
00089 !     Indices of the (method) grid cell in which the point was found.
00090 
00091       Logical, Intent (In)            :: msklocs (npoints)
00092 
00093 !     Mask of source locations array, only true srclocs need to be considered for target cells.
00094 
00095       Integer, Intent (In)            :: num_neigh
00096 !
00097 !     Last dimension of neighbors corner "neighbors_3d"
00098 !
00099       Real, Intent (In)               :: corner_x (corner_shape_3d (1,1,1):corner_shape_3d(2,1,1), 
00100                                                    corner_shape_3d (1,2,1):corner_shape_3d(2,2,1), 
00101                                                    corner_shape_3d (1,3,1):corner_shape_3d(2,3,1), 
00102                                                    nbr_corners(1))
00103 
00104       Real, Intent (In)               :: corner_y (corner_shape_3d (1,1,2):corner_shape_3d(2,1,2), 
00105                                                    corner_shape_3d (1,2,2):corner_shape_3d(2,2,2), 
00106                                                    corner_shape_3d (1,3,2):corner_shape_3d(2,3,2), 
00107                                                    nbr_corners(2))
00108 
00109       Real, Intent (In)               :: corner_z (corner_shape_3d (1,1,3):corner_shape_3d(2,1,3), 
00110                                                    corner_shape_3d (1,2,3):corner_shape_3d(2,2,3), 
00111                                                    corner_shape_3d (1,3,3):corner_shape_3d(2,3,3), 
00112                                                    nbr_corners(3))
00113 !
00114 !     Local corners
00115 !
00116 !
00117 ! !OUTPUT PARAMETERS:
00118 !
00119       Integer, Intent (Out)           :: nbr_cells(ncpl)
00120 !
00121 !     Number of cells found per target cell
00122 !
00123       Integer, Intent (Out)           :: ierror
00124 
00125 !     Returns the error code of PSMILE_Neigh_cells_irreg2_real;
00126 !             ierror = 0 : No error
00127 !             ierror > 0 : Severe error
00128 !
00129 ! !DEFINED PARAMETERS:
00130 !
00131 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00132 !               = 2 ** ndim_3d
00133 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00134 !               = 2 ** ndim_2d
00135 !
00136       Integer, Parameter              :: n_corners_3d = 2 * 2 * 2
00137       Integer, Parameter              :: n_corners_2d = 2 * 2
00138       Integer, Parameter              :: n_corners_1d = 2
00139 !
00140 ! !LOCAL VARIABLES
00141 !
00142 !     ... for local cell search
00143 !
00144       Real                           :: tgt_min, tgt_max
00145 
00146       Type (line_real), Allocatable  :: tgt(:)
00147       Type (line_real), Allocatable  :: src(:)
00148  
00149       Type (memo), Pointer :: stack(:)
00150       Type (memo), Pointer :: new_stack(:)
00151       Integer, Pointer     :: tmp_neighcells_3d(:,:)
00152       Integer, Pointer     :: new_neighcells_3d(:,:)
00153       Integer, Pointer     :: boundary_cell(:,:)
00154       Integer, Allocatable :: visited(:,:)
00155       Integer              :: unfolded_check
00156       Logical              :: overlap
00157 
00158       Integer              :: j1, j2
00159       Integer              :: ijdst(2,4)
00160       Integer              :: edge(2,4)
00161 
00162       Integer, Parameter   :: stacksize_inc = 512
00163       Integer, Parameter   :: bnd_size_inc  = 128
00164       Integer              :: nbr_cells_inc
00165 
00166       Integer              :: index
00167       Integer              :: stacksize
00168       Integer              :: bnd_size
00169 
00170       Integer                         :: i, j, n, m
00171       Integer                         :: n_corners
00172       Integer                         :: ii, jj, nn, mm, mold
00173 !
00174 !     ... For locations controlled
00175 !
00176       Integer                         :: length (ndim_3d)
00177 !
00178       Integer                         :: nca (ndim_3d)
00179       Integer                         :: nbr_cells_tot
00180 
00181       Integer                         :: grid_id
00182 
00183       Type(Grid), Pointer             :: gp
00184 !
00185 !     ... for error handling
00186 !
00187       Integer, Parameter              :: nerrp = 2
00188       Integer                         :: ierrp (nerrp)
00189 !
00190 #ifdef META_POST
00191       Integer                         :: rank
00192       Real                            :: shift_x, shift_y
00193       Integer, External               :: mp_fill_ov
00194       Integer, External               :: mp_fill_vi
00195       Integer, External               :: mp_draw_red
00196       Integer, External               :: mp_draw_black
00197       common / shift / shift_x, shift_y
00198 #endif
00199 !
00200 ! !DESCRIPTION:
00201 !
00202 ! Subroutine "PSMILe_Neigh_cells_3d_real" searches the cells on the donor
00203 ! grid that overlap with a given target cell.
00204 !
00205 ! Subroutine "PSMILe_Neigh_cells_3d_real" is designed for (target) grids
00206 ! of type "PRISM_irrlonlatvrt".
00207 !
00208 ! !TO DO:
00209 !
00210 !===> Check what we have to do for linear interpolation (interpolation_mode) in the vertical
00211 !      Is this already covered by a subsequent call to neigh_trili_3d?
00212 !
00213 !===> Cyclic boundary conditions
00214 !
00215 !
00216 ! !REVISION HISTORY:
00217 !
00218 !   Date      Programmer   Description
00219 ! ----------  ----------   -----------
00220 ! 06.06.19    R. Redler    created
00221 !
00222 !EOP
00223 !----------------------------------------------------------------------
00224 !
00225 !  $Id: psmile_neigh_cells_3d_real.F90 2082 2009-10-21 13:31:19Z hanke $
00226 !  $Author: hanke $
00227 !
00228    Character(len=len_cvs_string), save :: mycvs = 
00229        '$Id: psmile_neigh_cells_3d_real.F90 2082 2009-10-21 13:31:19Z hanke $'
00230 !
00231 !----------------------------------------------------------------------
00232 !
00233       data ijdst / -1, 0, 0, -1, 1, 0, 0, 1 /
00234 
00235       data edge / 1, 1, 2, 1, 2, 2, 1, 2 /
00236 !
00237 !  Relative indices of corner
00238 !
00239       data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00240 !
00241 !----------------------------------------------------------------------
00242 !
00243 !===> Prologue
00244 !
00245 #ifdef VERBOSE
00246       print 9990, trim(ch_id)
00247 
00248       call psmile_flushstd
00249 #endif /* VERBOSE */
00250 !
00251 !===> Control interpolation mode
00252 !
00253       if (interpolation_mode < ndim_2d .or. interpolation_mode > ndim_2d) then
00254          ierrp (1) = interpolation_mode
00255          ierror = PRISM_Error_Internal
00256 
00257          call psmile_error ( ierror, 'unsupported interpolation mode', &
00258                              ierrp, 1, __FILE__, __LINE__ )
00259          return
00260       endif
00261 !
00262 !===> Initialization
00263 !
00264       ierror = 0
00265       overlap = .false.
00266 
00267       shlon(0) =    0.0
00268       shlon(1) = -360.0
00269       shlon(2) =  360.0
00270 
00271       ! set an increment close to ncpl but an
00272       ! integer multiple of some power of two
00273       nbr_cells_inc = ((ncpl/16)+1)*16
00274 
00275       n_corners = nca (interpolation_mode)
00276 
00277       if ( n_corners /= N_CORNERS_2D ) then
00278          ierrp (1) = n_corners
00279          ierror = PRISM_Error_Internal
00280          call psmile_error ( ierror, 'unsupported number of corners', &
00281               ierrp, 1, __FILE__, __LINE__ )
00282       endif
00283 
00284       grid_id = use_how(3)
00285       gp => Grids(grid_id)
00286 
00287       length(1) = grid_shape(2,1) - grid_shape(1,1) + 1
00288       length(2) = grid_shape(2,2) - grid_shape(1,2) + 1
00289       length(3) = grid_shape(2,3) - grid_shape(1,3) + 1
00290       !
00291       !===> Find overlapping cells
00292       !
00293 #ifdef PRISM_ASSERTION
00294       do n = 1, ncpl
00295          if ( srclocs(1,n) < grid_shape(1,1) .or. &
00296               srclocs(1,n) > grid_shape(2,1) .or. &
00297               srclocs(2,n) < grid_shape(1,2) .or. &
00298               srclocs(2,n) > grid_shape(2,2) .or. &
00299               srclocs(3,n) < grid_shape(1,3) .or. &
00300               srclocs(3,n) > grid_shape(2,3)) exit
00301       enddo
00302 
00303       if (n <= ncpl) then
00304          print *, 'Incorrect index in n', n, srclocs (:, n)
00305          print *, 'grid shape ', grid_shape
00306          call psmile_assert (__FILE__, __LINE__, &
00307               'wrong index')
00308       endif
00309 #endif
00310 
00311 #ifdef META_POST
00312          rank = global_rank
00313 
00314          ! 595 points width for an A4 page, scaled by 180
00315          ! a shift my have to be applied in order to get
00316          ! the boxes on display
00317          if ( rank == META_RANK ) then
00318             open ( unit=99, file='cell_search.txt', form='formatted')
00319             write(99,'(a)') 'beginfig(-1)'
00320             write(99,'(a)') 'numeric u;'
00321             write(99,'(a3,f8.4,a1)') 'u:=', 595.0/180, ';'
00322             write(99,'(a)') 'pickup pencircle scaled .1pt;'
00323 
00324             shift_x = 0
00325             shift_y = 0
00326          endif
00327 #endif
00328       !
00329       !===>  Determine the i- and j-range on the donor grid for each target location
00330       !
00331       nbr_cells = 0
00332 
00333       m    = 0
00334       mm   = 0
00335       mold = 0
00336       nn   = 0
00337       
00338       allocate(visited(grid_shape(1,1):grid_shape(2,1), &
00339                        grid_shape(1,2):grid_shape(2,2)), STAT=ierror)
00340 
00341       if ( ierror > 0 ) then
00342          ierrp (1) = ierror
00343          ierrp (2) = length(1) * length(2)
00344          ierror = PRISM_Error_Alloc
00345          call psmile_error ( ierror, 'visited', &
00346               ierrp, 2, __FILE__, __LINE__ )
00347          return
00348       endif
00349 
00350       visited   = 0
00351 
00352       bnd_size = bnd_size_inc
00353 
00354       allocate(search%boundary_cell(bnd_size,5), STAT=ierror)
00355 
00356       if ( ierror > 0 ) then
00357          ierrp (1) = ierror
00358          ierrp (2) = 5*bnd_size
00359          ierror = PRISM_Error_Alloc
00360          call psmile_error ( ierror, 'boundary_cell', &
00361               ierrp, 2, __FILE__, __LINE__ )
00362          return
00363       endif
00364 
00365       search%boundary_cell(:,:) = PSMILe_Undef
00366 
00367       stacksize = stacksize_inc
00368 
00369       allocate(stack(stacksize), STAT=ierror)
00370 
00371       if ( ierror > 0 ) then
00372          ierrp (1) = ierror
00373          ierrp (2) = stacksize
00374          ierror = PRISM_Error_Alloc
00375          call psmile_error ( ierror, 'stack', &
00376               ierrp, 2, __FILE__, __LINE__ )
00377          return
00378       endif
00379 
00380       nbr_cells_tot = nbr_cells_inc
00381 
00382       allocate(tmp_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
00383 
00384       if ( ierror > 0 ) then
00385          ierrp (1) = ierror
00386          ierrp (2) = 2 * nbr_cells_tot
00387          ierror = PRISM_Error_Alloc
00388          call psmile_error ( ierror, 'tmp_neighcells_3d', &
00389               ierrp, 2, __FILE__, __LINE__ )
00390          return
00391       endif
00392 
00393       allocate(tgt(tgt_corners(1)), src(n_corners), STAT=ierror)
00394 
00395       if ( ierror > 0 ) then
00396          ierrp (1) = ierror
00397          ierrp (2) = 8
00398          ierror = PRISM_Error_Alloc
00399          call psmile_error ( ierror, 'tgt and src struct', &
00400               ierrp, 2, __FILE__, __LINE__ )
00401          return
00402       endif
00403 
00404       if ( gp%grid_type == PRISM_Gaussreduced_regvrt ) then
00405 
00406          do n = 1, npoints
00407 
00408             ! msklocs contains source cells that shall not be considered
00409             ! as starting cells for further local search. In particular
00410             ! this is required for reglonlatvrt source cell grids.
00411 
00412             if ( .not. msklocs(n) ) cycle
00413 
00414             m  = m + 1
00415 
00416             do j1 = 0, tgt_corners(1)-1
00417                j2 = mod(j1+1,tgt_corners(1))
00418                tgt(j1+1)%p1%x = tgt_cell(1)%vector(j1*ncpl+m)
00419                tgt(j1+1)%p1%y = tgt_cell(2)%vector(j1*ncpl+m)
00420                tgt(j1+1)%p2%x = tgt_cell(1)%vector(j2*ncpl+m)
00421                tgt(j1+1)%p2%y = tgt_cell(2)%vector(j2*ncpl+m)
00422             enddo
00423             !
00424             ! convert target longitudes to [0,360] interval if necessary
00425             !
00426             tgt_min =  999.0d0
00427             tgt_max = -999.0d0
00428 
00429             do j1 = 1, tgt_corners(1)
00430               tgt_min = min(tgt_min,tgt(j1)%p1%x)
00431               tgt_max = max(tgt_max,tgt(j1)%p1%x)
00432             enddo
00433 
00434 #ifdef DEBUGX
00435             !
00436             ! Find a particular entry
00437             !
00438             if ( tgt(1)%p1%y > -1 .and. tgt(1)%p1%y < 1 ) then
00439                print *, ' selected x Corners for ', m, tgt(:)%p1%x
00440                print *, ' selected y Corners for ', m, tgt(:)%p1%y
00441             endif
00442 #endif /* DEBUGX */
00443 
00444             if ( tgt_max - tgt_min > 180.0d0 ) then
00445 
00446                do j1 = 1, tgt_corners(1)
00447 
00448                   if ( tgt(j1)%p1%x > 360.0d0 ) then
00449                        tgt(j1)%p1%x = tgt(j1)%p1%x - 360.0d0
00450                        tgt_cell(1)%vector((j1-1)*ncpl+m) = &
00451                        tgt_cell(1)%vector((j1-1)*ncpl+m) - 360.0d0
00452                   endif
00453 
00454                   if ( tgt(j1)%p1%x < zero    ) then
00455                        tgt(j1)%p1%x = tgt(j1)%p1%x + 360.0d0
00456                        tgt_cell(1)%vector((j1-1)*ncpl+m) = &
00457                        tgt_cell(1)%vector((j1-1)*ncpl+m) + 360.0d0
00458                   endif
00459 
00460                   j2 = mod(j1,tgt_corners(1))
00461 
00462                   if ( tgt(j1)%p2%x > 360.0d0 ) then
00463                        tgt(j1)%p2%x = tgt(j1)%p2%x - 360.0d0
00464                        tgt_cell(1)%vector(j2*ncpl+m) = &
00465                        tgt_cell(1)%vector(j2*ncpl+m) - 360.0d0
00466                   endif
00467 
00468                   if ( tgt(j1)%p2%x < zero    ) then
00469                       tgt(j1)%p2%x = tgt(j1)%p2%x + 360.0d0
00470                       tgt_cell(1)%vector(j2*ncpl+m) = &
00471                       tgt_cell(1)%vector(j2*ncpl+m) + 360.0d0
00472                   endif
00473 
00474                enddo
00475 
00476             endif
00477             !
00478             ! For each new target cell the stack is reused and
00479             ! the first entry in the stack is initialised.
00480             !
00481             index = 1
00482 
00483             stack(index)%i   = srclocs(1,n)
00484             stack(index)%dir = 0
00485             !
00486             ! ... maybe we should first test whether our initial guess is really good
00487             !
00488             !
00489             ! Test current cell
00490             !
00491             ii = stack(index)%i
00492             !
00493             do j1 = 1, n_corners
00494                j2 = mod(j1,tgt_corners(1))+1
00495                src(j1)%p1%x = corner_x(ii, 1,1,edge(1,j1))
00496                src(j1)%p1%y = corner_y( 1,ii,1,edge(2,j1))
00497                src(j1)%p2%x = corner_x(ii, 1,1,edge(1,j2))
00498                src(j1)%p2%y = corner_y( 1,ii,1,edge(2,j2))
00499             enddo
00500 
00501             call psmile_overlap_real ( n_corners, tgt_corners(1), src, tgt, overlap )
00502 
00503             visited(ii,1) = m
00504 
00505             if ( overlap ) then
00506                !
00507                ! Store results of first cell
00508                !
00509                ! ... allocate new memory if necessary
00510                !
00511                if ( nn + 1 > nbr_cells_tot ) then
00512 
00513                   nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00514 
00515                   allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
00516 
00517                   if ( ierror > 0 ) then
00518                      ierrp (1) = ierror
00519                      ierrp (2) = 2 * nbr_cells_tot
00520                      ierror = PRISM_Error_Alloc
00521                      call psmile_error ( ierror, 'new_neighcells_3d', &
00522                           ierrp, 2, __FILE__, __LINE__ )
00523                      return
00524                   endif
00525 
00526                   new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
00527 
00528                   deallocate(tmp_neighcells_3d, STAT=ierror)
00529 
00530                   if ( ierror > 0 ) then
00531                      ierrp (1) = ierror
00532                      ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
00533                      ierror = PRISM_Error_Dealloc
00534                      call psmile_error ( ierror, 'tmp_neighcells_3d', &
00535                           ierrp, 2, __FILE__, __LINE__ )
00536                      return
00537                   endif
00538 
00539                   tmp_neighcells_3d => new_neighcells_3d
00540 
00541                endif
00542 
00543                nn = nn + 1
00544 
00545                tmp_neighcells_3d(1,nn) = stack(index)%i
00546                tmp_neighcells_3d(2,nn) = 1
00547                tmp_neighcells_3d(3,nn) = srclocs(3,n)
00548 
00549                nbr_cells(m) = 1
00550 
00551             endif
00552 
00553 #ifdef META_POST
00554 !rr         if ( mod(m,GRIDPOINT) == 0 ) then
00555             if ( m == GRIDPOINT .and. rank == META_RANK ) then
00556                
00557                shift_x = -tgt(1)%p1%x+180.0/6.0
00558                shift_y = -tgt(1)%p1%y+180.0/6.0
00559 
00560                ii = stack(index)%i
00561 
00562                if ( overlap ) then
00563                   ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00564 
00565                   do j1 = 1, n_corners
00566                      print *, ' source cell point o', ii, m, src(j1)%p1%x, src(j1)%p1%y
00567                      ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00568                   enddo
00569                   print *, ' --------- '
00570                else
00571                   ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00572 
00573                   do j1 = 1, n_corners
00574                      print *, ' source cell point v', ii, m, src(j1)%p1%x, src(j1)%p1%y
00575                      ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00576                   enddo
00577                   print *, ' --------- '
00578 
00579                endif
00580 
00581             endif
00582 #endif
00583             !
00584             ! Store further results with Andreas algorithm
00585             !
00586             do while ( index > 0 )
00587                !
00588                ! Change direction
00589                !
00590                ! dir = 1: east; dir = 2: south; dir = 3: west;, dir = 4: north
00591                !
00592                stack(index)%dir = stack(index)%dir + 1
00593 
00594                if ( stack(index)%dir > 4 ) then
00595                   !
00596                   ! All 4 directions are already tested
00597                   ! for this cell so we have to go back
00598                   !
00599                   index = index - 1
00600                   cycle
00601                endif
00602                !
00603                ! Go to the next cell (cyclic boundary conditions already in face)
00604                !
00605                ii = gp%face(stack(index)%i,2*stack(index)%dir)
00606                !
00607                ! Check that indices are in the allowed range
00608                !
00609                if ( ii < grid_shape(1,1) .or. ii > grid_shape(2,1) ) then
00610 
00611                   if ( m > mold ) then
00612 
00613                      mold = m
00614 
00615                      if ( mm + 1 > bnd_size ) then
00616 
00617                         bnd_size = bnd_size + bnd_size_inc
00618 
00619                         allocate(boundary_cell(bnd_size,5), STAT=ierror)
00620 
00621                         if ( ierror > 0 ) then
00622                            ierrp (1) = ierror
00623                            ierrp (2) = 5 * bnd_size
00624                            ierror = PRISM_Error_Alloc
00625                            call psmile_error ( ierror, 'boundary_cell', &
00626                                 ierrp, 2, __FILE__, __LINE__ )
00627                            return
00628                         endif
00629 
00630                         boundary_cell(1:mm,1:5) = search%boundary_cell(1:mm,1:5)
00631 
00632                         boundary_cell(mm+1:bnd_size,1:5) = PSMILe_Undef
00633 
00634                         deallocate(search%boundary_cell, STAT=ierror)
00635 
00636                         if ( ierror > 0 ) then
00637                            ierrp (1) = ierror
00638                            ierrp (2) = 5 * ( bnd_size - bnd_size_inc )
00639                            ierror = PRISM_Error_Dealloc
00640                            call psmile_error ( ierror, 'boundary_cell', &
00641                                 ierrp, 2, __FILE__, __LINE__ )
00642                            return
00643                         endif
00644 
00645                         search%boundary_cell => boundary_cell
00646 
00647                      endif
00648 
00649                      mm = mm + 1
00650 
00651                      search%boundary_cell(mm,1) = m
00652                      search%boundary_cell(mm,5) = psmile_rank
00653 
00654                   endif
00655 
00656                   cycle
00657                endif
00658                !
00659                ! When cell is already checked for this target ...
00660                !
00661                if ( visited(ii,1) == m ) cycle
00662                !
00663                ! Test current cell
00664                !
00665                do j1 = 1, n_corners
00666                   j2 = mod(j1,tgt_corners(1))+1
00667                   src(j1)%p1%x = corner_x(ii, 1,1,edge(1,j1))
00668                   src(j1)%p1%y = corner_y( 1,ii,1,edge(2,j1))
00669                   src(j1)%p2%x = corner_x(ii, 1,1,edge(1,j2))
00670                   src(j1)%p2%y = corner_y( 1,ii,1,edge(2,j2))
00671                enddo
00672 
00673                call psmile_overlap_real ( n_corners, tgt_corners(1), src, tgt, overlap )
00674                !
00675                ! Mark cell as checked for this cycle
00676                !
00677                visited(ii,1) = m
00678 #ifdef META_POST
00679 !rr            if ( mod(m,GRIDPOINT) == 0 ) then
00680                if ( m == GRIDPOINT .and. rank == META_RANK ) then
00681 
00682                   if ( overlap ) then
00683                      ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00684                      do j1 = 1,  n_corners
00685                         print *, ' source cell point o', ii, m, src(j1)%p1%x, src(j1)%p1%y
00686                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00687                      enddo
00688                   else
00689                      ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00690                      do j1 = 1,  n_corners
00691                         print *, ' source cell point v', ii, m, src(j1)%p1%x, src(j1)%p1%y
00692                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00693                      enddo
00694                   endif
00695                   print *, ' --------- '
00696 
00697                endif
00698 #endif
00699                if ( .not. overlap ) then
00700 
00701                   if ( stack(index)%dir == 1 .or. stack(index)%dir == 4 ) then
00702                      !
00703                      ! Make one more test to the west if not already visited as the
00704                      ! connectivity provided by the face array may have some gaps
00705                      !
00706                      jj = gp%face(ii,6)
00707                      !
00708                      if ( jj < grid_shape(1,1) .or. jj > grid_shape(2,1) ) cycle
00709                      !
00710                      if ( visited(jj,1)  == m ) jj = gp%face(ii,6)
00711                      !
00712                      if ( jj < grid_shape(1,1) .or. jj > grid_shape(2,1) ) cycle
00713                      !
00714                      if ( visited(jj,1) == m ) cycle
00715                      !
00716                      do j1 = 1, n_corners
00717                         j2 = mod(j1,tgt_corners(1))+1
00718                         src(j1)%p1%x = corner_x(jj, 1,1,edge(1,j1))
00719                         src(j1)%p1%y = corner_y( 1,jj,1,edge(2,j1))
00720                         src(j1)%p2%x = corner_x(jj, 1,1,edge(1,j2))
00721                         src(j1)%p2%y = corner_y( 1,jj,1,edge(2,j2))
00722                      enddo
00723 
00724                      call psmile_overlap_real ( n_corners, tgt_corners(1), src, tgt, overlap )
00725                      !
00726                      ! Mark cell as checked for this extra test
00727                      !
00728                      visited(jj,1) = m
00729 #ifdef META_POST
00730 !rr                  if ( mod(m,GRIDPOINT) == 0 ) then
00731                      if ( m == GRIDPOINT .and. rank == META_RANK ) then
00732 
00733                         if ( overlap ) then
00734                            ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00735                            do j1 = 1,  n_corners
00736                               print *, ' source cell point o', ii, m, src(j1)%p1%x, src(j1)%p1%y
00737                               ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00738                            enddo
00739                         else
00740                            ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00741                            do j1 = 1,  n_corners
00742                               print *, ' source cell point v', ii, m, src(j1)%p1%x, src(j1)%p1%y
00743                               ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00744                            enddo
00745                         endif
00746                         print *, ' --------- '
00747 
00748                      endif
00749 #endif
00750 
00751 #ifdef MAKE_A_TEST_TO_THE_EAST
00752                      if ( .not. overlap ) then
00753                         !
00754                         ! Make one more test to the east if not already visited as the
00755                         ! connectivity provided by the face array may have some gaps
00756                         !
00757                         jj = gp%face(ii,6)
00758                         !
00759                         if ( jj < grid_shape(1,1) .or. jj > grid_shape(2,1) ) cycle
00760                         !
00761                         if ( visited(jj,1) == m ) cycle
00762                         !
00763                         do j1 = 1, n_corners
00764                            j2 = mod(j1,tgt_corners(1))+1
00765                            src(j1)%p1%x = corner_x(jj, 1,1,edge(1,j1))
00766                            src(j1)%p1%y = corner_y( 1,jj,1,edge(2,j1))
00767                            src(j1)%p2%x = corner_x(jj, 1,1,edge(1,j2))
00768                            src(j1)%p2%y = corner_y( 1,jj,1,edge(2,j2))
00769                         enddo
00770 
00771                         call psmile_overlap_real ( n_corners, tgt_corners(1), src, tgt, overlap )
00772                         !
00773                         ! Mark cell as checked for this extra test
00774                         !
00775                         visited(jj,1) = m
00776 #ifdef META_POST
00777 !rr                     if ( mod(m,GRIDPOINT) == 0 ) then
00778                         if ( m == GRIDPOINT .and. rank == META_RANK ) then
00779 
00780                            if ( overlap ) then
00781                               ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00782                               do j1 = 1,  n_corners
00783                                  print *, ' source cell point o', ii, m, src(j1)%p1%x, src(j1)%p1%y
00784                                  ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00785                               enddo
00786                            else
00787                               ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00788                               do j1 = 1,  n_corners
00789                                  print *, ' source cell point v', ii, m, src(j1)%p1%x, src(j1)%p1%y
00790                                  ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00791                               enddo
00792                            endif
00793                            print *, ' --------- '
00794 
00795                         endif
00796 #endif
00797                      endif
00798 #endif
00799                      if ( overlap ) then
00800                         !
00801                         ! We found a new source cell
00802                         !
00803                         ! ... allocate new memory if necessary
00804                         !
00805                         if ( nn + 1 > nbr_cells_tot ) then
00806 
00807                            nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00808 
00809                            allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
00810 
00811                            if ( ierror > 0 ) then
00812                               ierrp (1) = ierror
00813                               ierrp (2) = 2 * nbr_cells_tot
00814                               ierror = PRISM_Error_Alloc
00815                               call psmile_error ( ierror, 'new_neighcells_3d', &
00816                                    ierrp, 2, __FILE__, __LINE__ )
00817                               return
00818                            endif
00819 
00820                            new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
00821 
00822                            deallocate(tmp_neighcells_3d, STAT=ierror)
00823 
00824                            if ( ierror > 0 ) then
00825                               ierrp (1) = ierror
00826                               ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
00827                               ierror = PRISM_Error_Dealloc
00828                               call psmile_error ( ierror, 'tmp_neighcells_3d', &
00829                                    ierrp, 2, __FILE__, __LINE__ )
00830                               return
00831                            endif
00832 
00833                            tmp_neighcells_3d => new_neighcells_3d
00834 
00835                         endif
00836                         !
00837                         ! Store results
00838                         !
00839                         nn = nn + 1
00840 
00841                         tmp_neighcells_3d(1,nn) = jj
00842                         tmp_neighcells_3d(2,nn) = 1
00843                         tmp_neighcells_3d(3,nn) = srclocs(3,n)
00844 
00845                         nbr_cells(m) = nbr_cells(m) + 1
00846 
00847                      endif
00848 
00849                   endif ! stack(index)%dir == 1 .or. stack(index)%dir == 4
00850 
00851                   cycle
00852 
00853                else
00854 
00855                   !
00856                   ! We found a new source cell
00857                   !
00858                   ! ... allocate new memory if necessary
00859                   !
00860                   if ( nn + 1 > nbr_cells_tot ) then
00861 
00862                      nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00863 
00864                      allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
00865 
00866                      if ( ierror > 0 ) then
00867                         ierrp (1) = ierror
00868                         ierrp (2) = 2 * nbr_cells_tot
00869                         ierror = PRISM_Error_Alloc
00870                         call psmile_error ( ierror, 'new_neighcells_3d', &
00871                              ierrp, 2, __FILE__, __LINE__ )
00872                         return
00873                      endif
00874 
00875                      new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
00876 
00877                      deallocate(tmp_neighcells_3d, STAT=ierror)
00878 
00879                      if ( ierror > 0 ) then
00880                         ierrp (1) = ierror
00881                         ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
00882                         ierror = PRISM_Error_Dealloc
00883                         call psmile_error ( ierror, 'tmp_neighcells_3d', &
00884                              ierrp, 2, __FILE__, __LINE__ )
00885                         return
00886                      endif
00887 
00888                      tmp_neighcells_3d => new_neighcells_3d
00889 
00890                   endif
00891                   !
00892                   ! Store results
00893                   !
00894                   nn = nn + 1
00895 
00896                   tmp_neighcells_3d(1,nn) = ii
00897                   tmp_neighcells_3d(2,nn) = 1
00898                   tmp_neighcells_3d(3,nn) = srclocs(3,n)
00899 
00900                   nbr_cells(m) = nbr_cells(m) + 1
00901 
00902                   !
00903                   ! Prepare for the next source cell to test
00904                   !
00905                   ! ... allocate new memory for the stack if necessary
00906                   !
00907                   if ( index + 1 > stacksize ) then
00908 
00909                      stacksize =  stacksize + stacksize_inc
00910 
00911                      allocate(new_stack(stacksize), STAT=ierror)
00912 
00913                      if ( ierror > 0 ) then
00914                         ierrp (1) = ierror
00915                         ierrp (2) = stacksize
00916                         ierror = PRISM_Error_Alloc
00917                         call psmile_error ( ierror, 'stack', &
00918                              ierrp, 2, __FILE__, __LINE__ )
00919                         return
00920                      endif
00921 
00922                      new_stack(1:index) = stack(1:index)
00923 
00924                      deallocate(stack, STAT=ierror)
00925 
00926                      if ( ierror > 0 ) then
00927                         ierrp (1) = ierror
00928                         ierrp (2) = stacksize - stacksize_inc
00929                         ierror = PRISM_Error_Dealloc
00930                         call psmile_error ( ierror, 'stack', &
00931                              ierrp, 2, __FILE__, __LINE__ )
00932                         return
00933                      endif
00934 
00935                      stack => new_stack
00936 
00937                   endif
00938                   !
00939                   ! ... initialise new stack entry
00940                   !
00941                   index = index + 1
00942                   stack(index)%i   = ii
00943                   stack(index)%dir = 0
00944 
00945                endif
00946 
00947             enddo
00948 #ifdef META_POST
00949 !rr         if ( mod(m,GRIDPOINT) == 0 ) then
00950             if ( m == GRIDPOINT .and. rank == META_RANK ) then
00951                do j1 = 1, tgt_corners(1)
00952                   ierror = mp_draw_black ( tgt(j1)%p1, tgt(j1)%p2 )
00953                   print *, ' target cell point ', m, tgt(j1)%p1%x, tgt(j1)%p1%y
00954                enddo
00955             endif
00956             ! for later graphical analysis print number of cells found for each target
00957             print *, ' nbr_cells(m) ', m, nbr_cells(m), n
00958 #endif
00959          enddo ! n loop over npoints
00960 
00961       else ! PRISM_Gaussreduced_regvrt
00962 
00963          !
00964          !-----------------------------------------------------
00965          ! branch for all Grids but Gaussreduced
00966          !-----------------------------------------------------
00967          !         
00968 
00969          do n = 1, npoints
00970 
00971             ! msklocs contains source cells that shall not be considered
00972             ! as starting cells for further local search. In particular
00973             ! this is required for reglonlatvrt source cell grids.
00974 
00975             if ( .not. msklocs(n) ) cycle
00976 
00977             m  = m + 1
00978 
00979             do j1 = 0, tgt_corners(1)-1
00980                j2 = mod(j1+1,tgt_corners(1))
00981                tgt(j1+1)%p1%x = tgt_cell(1)%vector(j1*ncpl+m)
00982                tgt(j1+1)%p1%y = tgt_cell(2)%vector(j1*ncpl+m)
00983                tgt(j1+1)%p2%x = tgt_cell(1)%vector(j2*ncpl+m)
00984                tgt(j1+1)%p2%y = tgt_cell(2)%vector(j2*ncpl+m)
00985             enddo
00986 #ifdef DEBUGX
00987             !
00988             ! Find a particular entry
00989             !
00990             if ( tgt(1)%p1%y > -1 .and. tgt(1)%p1%y < 1 ) then
00991                print *, ' selected x Corners for ', m, tgt(:)%p1%x
00992                print *, ' selected y Corners for ', m, tgt(:)%p1%y
00993             endif
00994 #endif /* DEBUGX */
00995             !
00996             ! convert target longitudes to [0,360] interval if necessary
00997             !
00998             tgt_min =  999.0d0
00999             tgt_max = -999.0d0
01000 
01001             do j1 = 1, tgt_corners(1)
01002               tgt_min = min(tgt_min,tgt(j1)%p1%x)
01003               tgt_max = max(tgt_max,tgt(j1)%p1%x)
01004             enddo
01005 
01006             if ( tgt_max - tgt_min > 180.0d0 ) then
01007 
01008                do j1 = 1, tgt_corners(1)
01009 
01010                   if ( tgt(j1)%p1%x > 360.0d0 ) then
01011                        tgt(j1)%p1%x = tgt(j1)%p1%x - 360.0d0
01012                        tgt_cell(1)%vector((j1-1)*ncpl+m) = &
01013                        tgt_cell(1)%vector((j1-1)*ncpl+m) - 360.0d0
01014                   endif
01015 
01016                   if ( tgt(j1)%p1%x < zero    ) then
01017                        tgt(j1)%p1%x = tgt(j1)%p1%x + 360.0d0
01018                       tgt_cell(1)%vector((j1-1)*ncpl+m) = &
01019                       tgt_cell(1)%vector((j1-1)*ncpl+m) + 360.0d0
01020                   endif
01021 
01022                   j2 = mod(j1,tgt_corners(1))
01023 
01024                   if ( tgt(j1)%p2%x > 360.0d0 ) then
01025                       tgt(j1)%p2%x = tgt(j1)%p2%x - 360.0d0
01026                       tgt_cell(1)%vector(j2*ncpl+m) = &
01027                       tgt_cell(1)%vector(j2*ncpl+m) - 360.0d0
01028                   endif
01029 
01030                   if ( tgt(j1)%p2%x < zero    ) then
01031                        tgt(j1)%p2%x = tgt(j1)%p2%x + 360.0d0
01032                        tgt_cell(1)%vector(j2*ncpl+m) = &
01033                        tgt_cell(1)%vector(j2*ncpl+m) + 360.0d0
01034                   endif
01035 
01036                enddo
01037 
01038             endif
01039 
01040             !
01041             ! For each new target cell the stack is reused and
01042             ! the first entry in the stack is initialised.
01043             !
01044             stack(:)%overlap = .false.
01045 
01046             index = 1
01047 
01048             stack(index)%i       = srclocs(1,n)
01049             stack(index)%j       = srclocs(2,n)
01050             stack(index)%dir     = 0
01051             stack(index)%overlap = .true.
01052             !
01053             ! ... maybe we should first test whether our initial guess is really good
01054             !
01055             visited(stack(index)%i,stack(index)%j) = m
01056             !
01057             ! Store results of first cell
01058             !
01059             ! ... allocate new memory if necessary
01060             !
01061             if ( nn + 1 > nbr_cells_tot ) then
01062 
01063                nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
01064 
01065                allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
01066 
01067                if ( ierror > 0 ) then
01068                   ierrp (1) = ierror
01069                   ierrp (2) = 2 * nbr_cells_tot
01070                   ierror = PRISM_Error_Alloc
01071                   call psmile_error ( ierror, 'new_neighcells_3d', &
01072                        ierrp, 2, __FILE__, __LINE__ )
01073                   return
01074                endif
01075 
01076                new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
01077 
01078                deallocate(tmp_neighcells_3d, STAT=ierror)
01079 
01080                if ( ierror > 0 ) then
01081                   ierrp (1) = ierror
01082                   ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
01083                   ierror = PRISM_Error_Dealloc
01084                   call psmile_error ( ierror, 'tmp_neighcells_3d', &
01085                        ierrp, 2, __FILE__, __LINE__ )
01086                   return
01087                endif
01088 
01089                tmp_neighcells_3d => new_neighcells_3d
01090 
01091             endif
01092 
01093             nn = nn + 1
01094 
01095             tmp_neighcells_3d(1,nn) = stack(index)%i
01096             tmp_neighcells_3d(2,nn) = stack(index)%j
01097             tmp_neighcells_3d(3,nn) = srclocs(3,n)
01098 
01099             nbr_cells(m) = 1
01100 
01101 #ifdef META_POST
01102 !rr         if ( mod(m,GRIDPOINT) == 0 ) then
01103             if ( m == GRIDPOINT .and. rank == META_RANK ) then
01104                
01105                shift_x = -tgt(1)%p1%x+180.0/6.0
01106                shift_y = -tgt(1)%p1%y+180.0/6.0
01107 
01108                ii = stack(index)%i
01109                jj = stack(index)%j
01110 
01111                if ( gp%grid_type == PRISM_Reglonlatvrt ) then
01112                   do j1 = 1, n_corners
01113                      j2 = mod(j1,tgt_corners(1))+1
01114                      src(j1)%p1%x = corner_x(ii, 1,1,edge(1,j1))
01115                      src(j1)%p1%y = corner_y( 1,jj,1,edge(2,j1))
01116                      src(j1)%p2%x = corner_x(ii, 1,1,edge(1,j2))
01117                      src(j1)%p2%y = corner_y( 1,jj,1,edge(2,j2))
01118                   enddo
01119                else
01120                   do j1 = 1, n_corners
01121                      j2 = mod(j1,tgt_corners(1))+1
01122                      src(j1)%p1%x = corner_x(ii,jj,1,j1)
01123                      src(j1)%p1%y = corner_y(ii,jj,1,j1)
01124                      src(j1)%p2%x = corner_x(ii,jj,1,j2)
01125                      src(j1)%p2%y = corner_y(ii,jj,1,j2)
01126                   enddo
01127                endif
01128 
01129                ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
01130 
01131                do j1 = 1, n_corners
01132                   print *, ' source cell point 1 o', ii, jj, nn+1, n, src(j1)%p1%x, src(j1)%p1%y
01133                   ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
01134                enddo
01135                print *, ' --------- '
01136             endif
01137 #endif
01138             !
01139             ! Store further results with Andreas algorithm
01140             !
01141 
01142             unfolded_check = 0
01143 
01144             do while ( index >= 0 )
01145 
01146                if ( index == 0 .and. unfolded_check > 0 ) then
01147 
01148                   call psmile_unfolded_check_dble( m,  & 
01149                                       unfolded_check,  &
01150                                       gp%grid_type,    &
01151                                       visited,         & 
01152                                       tgt_corners(1),  & 
01153                                       tgt_corners,     & 
01154                                       grid_shape,      & 
01155                                       corner_shape_3d, & 
01156                                       n_corners,       & 
01157                                       corner_x,        & 
01158                                       corner_y,        & 
01159                                       stacksize,       & 
01160                                       stack,           & 
01161                                       index,           & 
01162                                       ierror )
01163 
01164                   unfolded_check = 0
01165 
01166                   if ( index > 0 ) then
01167                      !
01168                      ! Store results of first cell
01169                      !
01170                      ! ... allocate new memory if necessary
01171                      !
01172                      if ( nn + 1 > nbr_cells_tot ) then
01173 
01174                         nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
01175 
01176                         allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
01177 
01178                         if ( ierror > 0 ) then
01179                            ierrp (1) = ierror
01180                            ierrp (2) = 2 * nbr_cells_tot
01181                            ierror = PRISM_Error_Alloc
01182                            call psmile_error ( ierror, 'new_neighcells_3d', &
01183                                 ierrp, 2, __FILE__, __LINE__ )
01184                            return
01185                         endif
01186 
01187                         new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
01188 
01189                         deallocate(tmp_neighcells_3d, STAT=ierror)
01190 
01191                         if ( ierror > 0 ) then
01192                            ierrp (1) = ierror
01193                            ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
01194                            ierror = PRISM_Error_Dealloc
01195                            call psmile_error ( ierror, 'tmp_neighcells_3d', &
01196                                 ierrp, 2, __FILE__, __LINE__ )
01197                            return
01198                         endif
01199 
01200                         tmp_neighcells_3d => new_neighcells_3d
01201 
01202                      endif
01203 
01204                      nn = nn + 1
01205 
01206                      tmp_neighcells_3d(1,nn) = stack(index)%i
01207                      tmp_neighcells_3d(2,nn) = stack(index)%j
01208                      tmp_neighcells_3d(3,nn) = srclocs(3,n)
01209 
01210                      nbr_cells(m) = nbr_cells(m) + 1
01211 
01212                      cycle
01213 
01214                   endif ! index > 0
01215 
01216                else if ( index == 0 ) then ! index == 0 .and. unfolded_check > 0
01217 
01218                   exit ! while loop
01219 
01220                endif ! index == 0 .and. unfolded_check > 0
01221 
01222                !
01223                ! Change direction
01224                !
01225                ! dir = 1: east; dir = 2: south; dir = 3: west;, dir = 4: north
01226                !
01227                stack(index)%dir = stack(index)%dir + 1
01228 
01229                if ( stack(index)%dir > 4 ) then
01230                   !
01231                   ! All 4 directions are already tested
01232                   ! for this cell so we have to go back
01233                   !
01234                   index = index - 1
01235                   cycle
01236                endif
01237                !
01238                ! Go to the next cell
01239                !
01240                ii = stack(index)%i + ijdst(1,stack(index)%dir)
01241                jj = stack(index)%j + ijdst(2,stack(index)%dir)
01242                !
01243                ! Apply cyclic boundary conditions if required
01244                !
01245                if ( cyclic(1) ) ii = mod(ii+length(1)-grid_shape(1,1),length(1))+grid_shape(1,1)
01246                if ( cyclic(2) ) jj = mod(jj+length(2)-grid_shape(1,2),length(2))+grid_shape(1,2)
01247                !
01248                ! Check that indices are in the allowed range
01249                !
01250                if ( ii < grid_shape(1,1) .or. ii > grid_shape(2,1) .or. &
01251                     jj < grid_shape(1,2) .or. jj > grid_shape(2,2) ) then
01252 
01253                   if ( gp%grid_type == PRISM_Irrlonlat_regvrt ) then
01254                      if ( jj < grid_shape(1,2) ) unfolded_check = 1
01255                      if ( jj > grid_shape(2,2) ) unfolded_check = 2
01256                   endif
01257 
01258                   if ( m > mold ) then
01259 
01260                      mold = m
01261 
01262                      if ( mm + 1 > bnd_size ) then
01263 
01264                         bnd_size = bnd_size + bnd_size_inc
01265 
01266                         allocate(boundary_cell(bnd_size,5), STAT=ierror)
01267 
01268                         if ( ierror > 0 ) then
01269                            ierrp (1) = ierror
01270                            ierrp (2) = 5 * bnd_size
01271                            ierror = PRISM_Error_Alloc
01272                            call psmile_error ( ierror, 'boundary_cell', &
01273                                 ierrp, 2, __FILE__, __LINE__ )
01274                            return
01275                         endif
01276 
01277                         boundary_cell(1:mm,1:5) = search%boundary_cell(1:mm,1:5)
01278 
01279                         boundary_cell(mm+1:bnd_size,1:5) = PSMILe_Undef
01280 
01281                         deallocate(search%boundary_cell, STAT=ierror)
01282 
01283                         if ( ierror > 0 ) then
01284                            ierrp (1) = ierror
01285                            ierrp (2) = 5 * ( bnd_size - bnd_size_inc )
01286                            ierror = PRISM_Error_Dealloc
01287                            call psmile_error ( ierror, 'boundary_cell', &
01288                                 ierrp, 2, __FILE__, __LINE__ )
01289                            return
01290                         endif
01291 
01292                         search%boundary_cell => boundary_cell
01293 
01294                      endif
01295 
01296                      mm = mm + 1
01297 
01298                      search%boundary_cell(mm,1) = m
01299                      search%boundary_cell(mm,5) = psmile_rank
01300 
01301                   endif
01302 
01303                   cycle
01304 
01305                endif
01306                !
01307                ! When cell is already checked for this target ...
01308                !
01309                if ( visited(ii,jj) == m ) cycle
01310                !
01311                ! Test current cell
01312                !
01313                if ( gp%grid_type == PRISM_Reglonlatvrt ) then
01314                   do j1 = 1, n_corners
01315                      j2 = mod(j1,tgt_corners(1))+1
01316                      src(j1)%p1%x = corner_x(ii, 1,1,edge(1,j1))
01317                      src(j1)%p1%y = corner_y( 1,jj,1,edge(2,j1))
01318                      src(j1)%p2%x = corner_x(ii, 1,1,edge(1,j2))
01319                      src(j1)%p2%y = corner_y( 1,jj,1,edge(2,j2))
01320                   enddo
01321                else
01322                   do j1 = 1, n_corners
01323                      j2 = mod(j1,tgt_corners(1))+1
01324                      src(j1)%p1%x = corner_x(ii,jj,1,j1)
01325                      src(j1)%p1%y = corner_y(ii,jj,1,j1)
01326                      src(j1)%p2%x = corner_x(ii,jj,1,j2)
01327                      src(j1)%p2%y = corner_y(ii,jj,1,j2)
01328                   enddo
01329                endif
01330 
01331                call psmile_overlap_real ( n_corners, tgt_corners(1), src, tgt, overlap )
01332                !
01333                ! Mark cell as checked for this cycle
01334                !
01335                visited(ii,jj) = m
01336 
01337                stack(index)%overlap = overlap
01338 #ifdef META_POST
01339 !rr            if ( mod(m,GRIDPOINT) == 0 ) then
01340                if ( m == GRIDPOINT .and. rank == META_RANK ) then
01341 
01342                   if ( overlap ) then
01343                      ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
01344                      do j1 = 1,  n_corners
01345                         print *, ' source cell point 2 o', ii, jj, nn+1, n, src(j1)%p1%x, src(j1)%p1%y
01346                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
01347                      enddo
01348                   else
01349                      ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
01350                      do j1 = 1,  n_corners
01351                         print *, ' source cell point v', ii, jj, nn+1, n, src(j1)%p1%x, src(j1)%p1%y
01352                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
01353                      enddo
01354                   endif
01355                   print *, ' --------- '
01356 
01357                endif
01358 #endif
01359 
01360                if ( .not. overlap ) then
01361 
01362                   cycle
01363 
01364                else
01365                   !
01366                   ! We found a new source cell
01367                   !
01368                   ! ... allocate new memory if necessary
01369                   !
01370                   if ( nn + 1 > nbr_cells_tot ) then
01371 
01372                      nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
01373 
01374                      allocate(new_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
01375 
01376                      if ( ierror > 0 ) then
01377                         ierrp (1) = ierror
01378                         ierrp (2) = 2 * nbr_cells_tot
01379                         ierror = PRISM_Error_Alloc
01380                         call psmile_error ( ierror, 'new_neighcells_3d', &
01381                              ierrp, 2, __FILE__, __LINE__ )
01382                         return
01383                      endif
01384 
01385                      new_neighcells_3d(1:3,1:nn) = tmp_neighcells_3d(1:3,1:nn)
01386 
01387                      deallocate(tmp_neighcells_3d, STAT=ierror)
01388 
01389                      if ( ierror > 0 ) then
01390                         ierrp (1) = ierror
01391                         ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
01392                         ierror = PRISM_Error_Dealloc
01393                         call psmile_error ( ierror, 'tmp_neighcells_3d', &
01394                              ierrp, 2, __FILE__, __LINE__ )
01395                         return
01396                      endif
01397 
01398                      tmp_neighcells_3d => new_neighcells_3d
01399 
01400                   endif
01401                   !
01402                   ! Store results
01403                   !
01404                   nn = nn + 1
01405 
01406                   tmp_neighcells_3d(1,nn) = ii ! - grid_shape(1,1) + 1
01407                   tmp_neighcells_3d(2,nn) = jj ! - grid_shape(1,2) + 1
01408                   tmp_neighcells_3d(3,nn) = srclocs(3,n) ! - grid_shape(1,3) + 1
01409 
01410                   nbr_cells(m) = nbr_cells(m) + 1
01411 
01412                   !
01413                   ! Prepare for the next source cell to test
01414                   !
01415                   ! ... allocate new memory for the stack if necessary
01416                   !
01417                   if ( index + 1 > stacksize ) then
01418 
01419                      stacksize =  stacksize + stacksize_inc
01420 
01421                      allocate(new_stack(stacksize), STAT=ierror)
01422 
01423                      if ( ierror > 0 ) then
01424                         ierrp (1) = ierror
01425                         ierrp (2) = stacksize
01426                         ierror = PRISM_Error_Alloc
01427                         call psmile_error ( ierror, 'stack', &
01428                              ierrp, 2, __FILE__, __LINE__ )
01429                         return
01430                      endif
01431 
01432                      new_stack(:)%overlap = .false.
01433 
01434                      new_stack(1:index) = stack(1:index)
01435 
01436                      deallocate(stack, STAT=ierror)
01437 
01438                      if ( ierror > 0 ) then
01439                         ierrp (1) = ierror
01440                         ierrp (2) = stacksize - stacksize_inc
01441                         ierror = PRISM_Error_Dealloc
01442                         call psmile_error ( ierror, 'stack', &
01443                              ierrp, 2, __FILE__, __LINE__ )
01444                         return
01445                      endif
01446 
01447                      stack => new_stack
01448 
01449                   endif
01450                   !
01451                   ! ... initialise new stack entry
01452                   !
01453                   index = index + 1
01454                   stack(index)%i   = ii
01455                   stack(index)%j   = jj
01456                   stack(index)%dir = 0
01457 
01458                endif
01459 
01460             enddo ! while
01461 
01462 #ifdef META_POST
01463 !rr         if ( mod(m,GRIDPOINT) == 0 ) then
01464             if ( m == GRIDPOINT .and. rank == META_RANK ) then
01465                do j1 = 1, tgt_corners(1)
01466                   ierror = mp_draw_black ( tgt(j1)%p1, tgt(j1)%p2 )
01467                   print *, ' target cell point ', m, tgt(j1)%p1%x, tgt(j1)%p1%y
01468                enddo
01469             endif
01470             ! for later graphical analysis print number of cells found for each target
01471             print *, ' nbr_cells(m) ', m, nbr_cells(m), n
01472 #endif
01473          enddo ! n loop over npoints
01474 
01475       endif ! Gaussreduced_regvrt
01476 
01477 #ifdef META_POST
01478       if ( rank == META_RANK ) then
01479          write(99,'(a)') 'endfig'
01480          write(99,'(a)') 'end'
01481          close(unit=99)
01482       endif
01483 #endif
01484 
01485       nbr_cells_tot = nn
01486       !
01487       ! Copy temporary results to final destination and deallocate temporary buffer
01488       !
01489       select case ( gp%grid_type )
01490 
01491       case ( PRISM_Reglonlatvrt )
01492 
01493          Allocate ( neighcells_3d(ndim_3d,nbr_cells_tot,n_corners), STAT=ierror)
01494 
01495          if (ierror /= 0) then
01496             ierrp (1) = 3 * nbr_cells_tot * n_corners
01497             ierror = PRISM_Error_Alloc
01498             call psmile_error ( ierror, 'neighcells_3d', &
01499                  ierrp, 1, __FILE__, __LINE__ )
01500             return
01501          endif
01502 
01503          neighcells_3d(1:3,1:nbr_cells_tot,1) = tmp_neighcells_3d(1:3,1:nbr_cells_tot)
01504 
01505          neighcells_3d(1,1:nbr_cells_tot,2) = neighcells_3d(1,1:nbr_cells_tot,1) + 1
01506          neighcells_3d(1,1:nbr_cells_tot,3) = neighcells_3d(1,1:nbr_cells_tot,1) + 1
01507          neighcells_3d(1,1:nbr_cells_tot,4) = neighcells_3d(1,1:nbr_cells_tot,1)
01508 
01509          neighcells_3d(2,1:nbr_cells_tot,2) = neighcells_3d(2,1:nbr_cells_tot,1)
01510          neighcells_3d(2,1:nbr_cells_tot,3) = neighcells_3d(2,1:nbr_cells_tot,1) + 1
01511          neighcells_3d(2,1:nbr_cells_tot,4) = neighcells_3d(2,1:nbr_cells_tot,1) + 1
01512 
01513          neighcells_3d(3,1:nbr_cells_tot,2) = neighcells_3d(3,1:nbr_cells_tot,1)
01514          neighcells_3d(3,1:nbr_cells_tot,3) = neighcells_3d(3,1:nbr_cells_tot,1)
01515          neighcells_3d(3,1:nbr_cells_tot,4) = neighcells_3d(3,1:nbr_cells_tot,1)
01516          !
01517          !===> Store other corner indices for corners 5, 6, 7, and 8 when required, e.g.
01518          !     for 2D1D interpolation with 1D in the vertical.
01519          !
01520          if ( n_corners == 8 ) then
01521             do n = 5, n_corners
01522                neighcells_3d(1,1:nbr_cells_tot,n) = neighcells_3d(1,1:nbr_cells_tot,n-4)
01523                neighcells_3d(2,1:nbr_cells_tot,n) = neighcells_3d(2,1:nbr_cells_tot,n-4)
01524                neighcells_3d(3,1:nbr_cells_tot,n) = neighcells_3d(3,1:nbr_cells_tot,1)+1
01525             enddo
01526          endif
01527          !
01528          !===> Set indices for cyclic coordinate directions
01529          !
01530          do i = 1, ndim_3d
01531 
01532             if (cyclic(i) .and. length(i) > 1) then
01533 
01534                do n = 1, n_corners
01535 
01536                   do j = 1, nbr_cells_tot
01537 
01538                      if ( neighcells_3d(i,j,n) < grid_shape(1,i)) then
01539                           neighcells_3d(i,j,n) = neighcells_3d(i,j,n) + length (i)
01540                      else if ( neighcells_3d(i,j,n) > grid_shape(2,i)) then
01541                           neighcells_3d(i,j,n) = neighcells_3d(i,j,n) - length (i)
01542                      endif
01543 
01544                   enddo
01545 
01546                end do ! n
01547 
01548             endif
01549 
01550          end do ! j
01551          !
01552          !===> Control special cases (2d hyperplanes)
01553          !
01554          do i = 1, ndim_3d
01555             if (length(i) == 1) then
01556                neighcells_3d(i,1:nbr_cells_tot,:) = grid_shape(1,i)
01557             endif
01558          end do
01559 
01560       case DEFAULT
01561 
01562          allocate(neighcells_3d(3,nbr_cells_tot,1), STAT=ierror)
01563 
01564          if ( ierror > 0 ) then
01565             ierrp (1) = ierror
01566             ierrp (2) = 3 * nbr_cells_tot
01567             ierror = PRISM_Error_Alloc
01568             call psmile_error ( ierror, 'neighcells_3d', &
01569                  ierrp, 2, __FILE__, __LINE__ )
01570             return
01571          endif
01572 
01573          neighcells_3d(1:3,1:nn,1) = tmp_neighcells_3d(1:3,1:nn)
01574 
01575       end select
01576 
01577       bnd_size = mm
01578 
01579       if ( bnd_size > 0 ) then
01580 
01581          allocate(boundary_cell(bnd_size,5), STAT=ierror)
01582 
01583          if ( ierror > 0 ) then
01584             ierrp (1) = ierror
01585             ierrp (2) = 5 * bnd_size
01586             ierror = PRISM_Error_Alloc
01587             call psmile_error ( ierror, 'boundary_cell', &
01588                  ierrp, 2, __FILE__, __LINE__ )
01589             return
01590          endif
01591 
01592          boundary_cell(1:bnd_size,1:5) = search%boundary_cell(1:bnd_size,1:5)
01593 
01594          ! We should memorise a cell only when neighbours have been found
01595 
01596          do i = 1, bnd_size
01597             if ( search%boundary_cell(i,1) /= PSMILe_Undef ) then
01598                if (nbr_cells(search%boundary_cell(mm,1)) == 0 ) &
01599                     boundary_cell(i,1) = PSMILe_Undef
01600             endif
01601 #ifdef DEBUGX
01602             !
01603             ! boundary cell indices 2, 3 and 4 are set in psmile_info_trs_loc_*
01604             !
01605             if ( search%boundary_cell(i,1) /= PSMILe_Undef ) then
01606                 print '(a,4i8)', 'boundary cell info ', &
01607                    i, boundary_cell(i,1), boundary_cell(i,5), nbr_cells(boundary_cell(i,1))
01608             endif
01609 #endif /* DEBUGX */
01610          enddo
01611 
01612          deallocate(search%boundary_cell, STAT=ierror)
01613 
01614          if ( ierror > 0 ) then
01615             ierrp (1) = ierror
01616             ierrp (2) = 5 * ( bnd_size - bnd_size_inc )
01617             ierror = PRISM_Error_Dealloc
01618             call psmile_error ( ierror, 'boundary_cell', &
01619                  ierrp, 2, __FILE__, __LINE__ )
01620             return
01621          endif
01622 
01623          search%boundary_cell => boundary_cell
01624 
01625       else
01626 
01627          deallocate(search%boundary_cell, STAT=ierror)
01628 
01629          if ( ierror > 0 ) then
01630             ierrp (1) = ierror
01631             ierrp (2) = 3 * nbr_cells_tot
01632             ierror = PRISM_Error_Dealloc
01633             call psmile_error ( ierror, 'search%boundary_cell', &
01634                  ierrp, 2, __FILE__, __LINE__ )
01635             return
01636          endif
01637 
01638       endif
01639 !
01640 !===> Deallocate local memory
01641 !
01642       deallocate(tgt, src, STAT=ierror)
01643 
01644       if ( ierror > 0 ) then
01645          ierrp (1) = ierror
01646          ierrp (2) = 8
01647          ierror = PRISM_Error_Dealloc
01648          call psmile_error ( ierror, 'tgt and src struct', &
01649               ierrp, 2, __FILE__, __LINE__ )
01650          return
01651       endif
01652 
01653       deallocate(tmp_neighcells_3d, STAT=ierror)
01654 
01655       if ( ierror > 0 ) then
01656          ierrp (1) = ierror
01657          ierrp (2) = 3 * nbr_cells_tot
01658          ierror = PRISM_Error_Dealloc
01659          call psmile_error ( ierror, 'tmp_neighcells_3d', &
01660               ierrp, 2, __FILE__, __LINE__ )
01661          return
01662       endif
01663 
01664       deallocate(stack, visited, STAT=ierror)
01665 
01666       if ( ierror > 0 ) then
01667          ierrp (1) = ierror
01668          ierrp (2) = stacksize
01669          ierror = PRISM_Error_Dealloc
01670          call psmile_error ( ierror, 'stack and visited', &
01671               ierrp, 2, __FILE__, __LINE__ )
01672          return
01673       endif
01674 !
01675 !===> All done
01676 !
01677 
01678 #ifdef VERBOSE
01679       print 9980, trim(ch_id)
01680 
01681       call psmile_flushstd
01682 #endif /* VERBOSE */
01683 !
01684 !  Formats:
01685 !
01686 9990 format (1x, a, ': psmile_neigh_cells_3d_real')
01687 9980 format (1x, a, ': psmile_neigh_cells_3d_real: eof')
01688 
01689       end subroutine psmile_neigh_cells_3d_real
01690 
01691       
01692 #ifdef META_POST
01693       integer function mp_fill_ov ( s1, s2, s3, s4 )
01694          use PSMILe, only : point_real
01695          Type(point_real) :: p1, p2, p3, p4
01696          Type(point_real) :: s1, s2, s3, s4
01697          Real             :: shift_x, shift_y
01698          common / shift / shift_x, shift_y
01699 
01700          p1 = s1
01701          p2 = s2
01702          p3 = s3
01703          p4 = s4
01704 
01705          p1%x = s1%x + shift_x
01706          p2%x = s2%x + shift_x
01707          p3%x = s3%x + shift_x
01708          p4%x = s4%x + shift_x
01709 
01710          p1%y = s1%y + shift_y
01711          p2%y = s2%y + shift_y
01712          p3%y = s3%y + shift_y
01713          p4%y = s4%y + shift_y
01714 
01715          write(99,'(a6,4(f9.3,a2,f9.3,a5),a27)') 'fill (',p1%x,'u,',p1%y,'u)--(' &
01716                                                          ,p2%x,'u,',p2%y,'u)--(' &
01717                                                          ,p3%x,'u,',p3%y,'u)--(' &
01718                                                          ,p4%x,'u,',p4%y,'u)','--cycle withcolor .6*white;'
01719          mp_fill_ov = 1
01720 
01721       end function
01722 
01723       integer function mp_fill_vi ( s1, s2, s3, s4 )
01724          use PSMILe, only : point_real
01725          Type(point_real) :: p1, p2, p3, p4
01726          Type(point_real) :: s1, s2, s3, s4
01727          Real             :: shift_x, shift_y
01728          common / shift / shift_x, shift_y
01729 
01730          p1 = s1
01731          p2 = s2
01732          p3 = s3
01733          p4 = s4
01734 
01735          p1%x = s1%x + shift_x
01736          p2%x = s2%x + shift_x
01737          p3%x = s3%x + shift_x
01738          p4%x = s4%x + shift_x
01739 
01740          p1%y = s1%y + shift_y
01741          p2%y = s2%y + shift_y
01742          p3%y = s3%y + shift_y
01743          p4%y = s4%y + shift_y
01744 
01745          write(99,'(a6,4(f9.3,a2,f9.3,a5),a27)') 'fill (',p1%x,'u,',p1%y,'u)--(' &
01746                                                          ,p2%x,'u,',p2%y,'u)--(' &
01747                                                          ,p3%x,'u,',p3%y,'u)--(' &
01748                                                          ,p4%x,'u,',p4%y,'u)','--cycle withcolor .9*white;'
01749          mp_fill_vi = 1
01750 
01751       end function
01752 
01753       integer function mp_draw_red ( s1, s2 )
01754          use PSMILe, only : point_real
01755          Type(point_real) :: s1, s2
01756          Type(point_real) :: p1, p2
01757          Real             :: shift_x, shift_y
01758          common / shift / shift_x, shift_y
01759 
01760          p1 = s1
01761          p2 = s2
01762          p1%x = s1%x + shift_x
01763          p2%x = s2%x + shift_x
01764 
01765          p1%y = s1%y + shift_y
01766          p2%y = s2%y + shift_y
01767 
01768          write(99,'(a6,2(f9.3,a2,f9.3,a5),a25)') 'draw (',p1%x,'u,',p1%y,'u)--(' &
01769                                                          ,p2%x,'u,',p2%y,'u)','--cycle withcolor red;'
01770          mp_draw_red = 1
01771 
01772       end function
01773 
01774       integer function mp_draw_black ( s1, s2 )
01775          use PSMILe, only : point_real
01776          Type(point_real) :: s1, s2
01777          Type(point_real) :: p1, p2
01778          Real             :: shift_x, shift_y
01779          common / shift / shift_x, shift_y
01780 
01781          p1 = s1
01782          p2 = s2
01783 
01784          p1%x = s1%x + shift_x
01785          p2%x = s2%x + shift_x
01786 
01787          p1%y = s1%y + shift_y
01788          p2%y = s2%y + shift_y
01789 
01790          write(99,'(a6,2(f9.3,a2,f9.3,a5),a27)') 'draw (',p1%x,'u,',p1%y,'u)--(' &
01791                                                          ,p2%x,'u,',p2%y,'u)','--cycle withcolor black;'
01792          mp_draw_black = 1
01793 
01794       end function
01795 #endif
Generated on 18 Mar 2011 for Oasis4 by  doxygen 1.6.1