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