psmile_neigh_cells_irreg2_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 #define GRIDPOINT 2345
00008 #undef  META_POST
00009 !
00010 ! !ROUTINE: PSMILe_Neigh_cells_irreg2_dble
00011 !
00012 ! !INTERFACE:
00013 
00014       subroutine psmile_neigh_cells_irreg2_dble ( use_how,      &
00015                grid_shape, interpolation_mode, cyclic,          &
00016                corner_shape_3d, nbr_corners,                    &
00017                corner_x, corner_y,                              &
00018                search,  control, tgt_cell, tgt_corners,         &
00019                npoints, srclocs, msklocs, ncpl,                 &
00020                num_neigh, nbr_cells, ierror )
00021 !
00022 ! !USES:
00023 !
00024       use PRISM_constants
00025 !
00026       use PSMILe, dummy_interface => PSMILe_Neigh_cells_irreg2_dble
00027 
00028       implicit none
00029 !
00030 ! !INPUT PARAMETERS:
00031 !
00032       Integer, Intent (In)            :: use_how(3)
00033 !
00034 !     container of information  under consideration
00035 !
00036       Integer, Intent (In)            :: grid_shape (2, ndim_3d)
00037 !
00038 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00039 !
00040       Integer, Intent (In)            :: interpolation_mode
00041 !
00042 !     Specifies the interpolation mode with
00043 !        interpolation_mode = 3 : 3D conservative remapping
00044 !                                 (not supported)
00045 !        interpolation_mode = 2 : 2D conservative remapping
00046 !        interpolation_mode = 1 : not supported
00047 !                                 (determines 2 corners)
00048 !
00049       Logical, Intent (In)            :: cyclic (ndim_3d)
00050 !    
00051 !     Does the block have cyclic coordinates in the corresponding direction ?
00052 !
00053       Integer, Intent (In)            :: corner_shape_3d(2,ndim_3d,ndim_3d)
00054 !
00055 !     Shape of local corners -> corner_x, corner_y
00056 !
00057       Integer, Intent (In)            :: nbr_corners(ndim_3d)
00058 !
00059 !     Number of corners of local corners corner_x, corner_y
00060 !
00061       Type (Enddef_search), Intent (InOut) :: search
00062 !
00063 !     Info's on coordinates to be searched
00064 !
00065       Integer, Intent (In)            :: control (2, ndim_3d, search%npart)
00066 !
00067 !     Control range
00068 !
00069       Type (dble_vector), Intent (InOut) :: tgt_cell(ndim_3d)
00070 !
00071 !     Cells for which we search
00072 !
00073       Integer, Intent (In)            :: tgt_corners
00074 !
00075 !     Number of corner in the target cell
00076 !
00077       Integer, Intent (In)            :: npoints (ndim_2d, search%npart)
00078 !
00079 !     Total number of locations to be transferred
00080 !
00081       Integer, Intent (In)            :: ncpl
00082 !
00083 !     Total number of locations to be transferred
00084 !
00085       Type (integer_vector), Intent (In) :: srclocs (2,search%npart)
00086 !
00087 !     Indices of the (method) grid cell in which the point was found.
00088 
00089       Type (logical_vector), Intent (In) :: msklocs (2,search%npart)
00090 
00091 !     Mask of source locations array, only true srclocs need to be considered for target cells.
00092 
00093       Integer, Intent (In)            :: num_neigh
00094 !
00095 !     Last dimension of neighbors corner "neighcells_3d"
00096 !
00097       Double Precision, Intent (In)   :: corner_x (corner_shape_3d (1,1,1):corner_shape_3d(2,1,1), 
00098                                                    corner_shape_3d (1,2,1):corner_shape_3d(2,2,1), 
00099                                                    nbr_corners(1))
00100 
00101       Double Precision, Intent (In)   :: corner_y (corner_shape_3d (1,1,2):corner_shape_3d(2,1,2), 
00102                                                    corner_shape_3d (1,2,2):corner_shape_3d(2,2,2), 
00103                                                    nbr_corners(2))
00104 !
00105 !     Local corners
00106 !
00107 !
00108 ! !OUTPUT PARAMETERS:
00109 !
00110       Integer, Intent (Out)           :: nbr_cells(ncpl)
00111 !
00112 !     Number of cells found per target cell
00113 !
00114       Integer, Intent (Out)           :: ierror
00115 
00116 !     Returns the error code of PSMILE_Neigh_cells_irreg2_dble;
00117 !             ierror = 0 : No error
00118 !             ierror > 0 : Severe error
00119 !
00120 ! !DEFINED PARAMETERS:
00121 !
00122 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00123 !               = 2 ** ndim_3d
00124 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00125 !               = 2 ** ndim_2d
00126 !
00127       Integer, parameter              :: n_corners_3d = 2 * 2 * 2
00128       Integer, parameter              :: n_corners_2d = 2 * 2
00129       Integer, parameter              :: n_corners_1d = 2
00130 !
00131 ! !LOCAL VARIABLES
00132 !
00133 !     ... for corner locations searched
00134 !
00135       Double Precision                :: tgt_min, tgt_max
00136 
00137       Type (line_dble), Allocatable   :: tgt(:)
00138       Type (line_dble), Allocatable   :: src(:)
00139  
00140       Type (memo), Pointer            :: stack(:)
00141       Type (memo), Pointer            :: new_stack(:)
00142       Integer, Pointer                :: tmp_neighcells_3d(:,:)
00143       Integer, Pointer                :: new_neighcells_3d(:,:)
00144       Integer, Pointer                :: boundary_cell(:,:)
00145       Integer, Allocatable            :: visited(:,:)
00146       Integer                         :: unfolded_check
00147       Logical                         :: overlap
00148 
00149       Integer                         :: j1, j2
00150       Integer                         :: ijdst(2,4)
00151       Integer                         :: edge(2,4)
00152 
00153       Integer, Parameter              :: stacksize_inc = 512
00154       Integer, Parameter              :: bnd_size_inc  = 128
00155       Integer                         :: nbr_cells_inc
00156 
00157       Integer                         :: index
00158       Integer                         :: stacksize
00159       Integer                         :: bnd_size
00160 
00161       Integer                         :: ipart
00162       Integer                         :: i, k, n, m, mm, mold
00163       Integer                         :: nend, nprev
00164       Integer                         :: lenij, lenijk, n_corners
00165       Integer                         :: ii, jj, nn, offset1, offset2
00166 !
00167 !     ... For locations controlled
00168 !
00169       Integer                         :: length (ndim_3d)
00170 !
00171       Integer                         :: ntgtlocs (ndim_2d)
00172       Integer                         :: nca (ndim_3d)
00173       Integer                         :: ndim ! nbr of dims for conservative remapping
00174       Integer                         :: nbr_cells_tot
00175 
00176       Integer                         :: grid_id
00177 
00178       Type(Grid), Pointer             :: gp
00179 !
00180 !     ... for error handling
00181 !
00182       Integer, parameter              :: nerrp = 2
00183       Integer                         :: ierrp (nerrp)
00184 !
00185 #ifdef META_POST
00186       Double Precision                :: shift_x, shift_y
00187       Integer, External               :: mp_fill_ov
00188       Integer, External               :: mp_fill_vi
00189       Integer, External               :: mp_draw_red
00190       Integer, External               :: mp_draw_black
00191       common / shift / shift_x, shift_y
00192 #endif
00193 !
00194 ! !DESCRIPTION:
00195 !
00196 ! Subroutine "PSMILe_Neigh_cells_irreg2_dble" searches the cells
00197 ! for the selected interpolation type (conserve2D respectively)
00198 ! on the the source-grid cells for the subgrid corners sent by
00199 ! the destination process.
00200 !
00201 ! Subroutine "PSMILe_Neigh_cells_irreg2_dble" is designed for
00202 ! (target) grids of type "PRISM_Irrlonlat_regvrt".
00203 !
00204 ! !REVISION HISTORY:
00205 !
00206 !   Date      Programmer   Description
00207 ! ----------  ----------   -----------
00208 ! 06.09.06    R. Redler    created
00209 !
00210 !EOP
00211 !----------------------------------------------------------------------
00212 !
00213 !  $Id: psmile_neigh_cells_irreg2_dble.F90 2082 2009-10-21 13:31:19Z hanke $
00214 !  $Author: hanke $
00215 !
00216    Character(len=len_cvs_string), save :: mycvs = 
00217        '$Id: psmile_neigh_cells_irreg2_dble.F90 2082 2009-10-21 13:31:19Z hanke $'
00218 !
00219 !----------------------------------------------------------------------
00220 !
00221      data ijdst / -1, 0, 0, -1, 1, 0, 0, 1 /
00222 
00223      data edge / 1, 1, 2, 1, 2, 2, 1, 2 /
00224 !
00225 !  Relative indices of corner
00226 !
00227      data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00228 !
00229 !----------------------------------------------------------------------
00230 !
00231 !  Initialization
00232 !
00233 #ifdef VERBOSE
00234       print 9990, trim(ch_id)
00235 
00236       call psmile_flushstd
00237 #endif /* VERBOSE */
00238 
00239 !
00240 !===> Control interpolation mode
00241 !
00242       if (interpolation_mode < 1 .or. interpolation_mode > ndim_3d) then
00243          ierrp (1) = interpolation_mode
00244          ierror = PRISM_Error_Internal
00245 
00246          call psmile_error ( ierror, &
00247             'unsupported interpolation mode in psmile_neigh_cells_irreg2_dble', &
00248             ierrp, 1, __FILE__, __LINE__ )
00249          return
00250       endif
00251 !
00252 !===> Initialization
00253 !
00254 !  lenijk = Number of entries in all   3 directions
00255 !
00256       ierror = 0
00257 
00258       ! set an increment close to ncpl but an
00259       ! integer multiple of some power of two
00260       nbr_cells_inc = max(16,floor(real(ncpl)/real(16))*16)
00261       nbr_cells_inc = ((ncpl/16)+1)*16
00262 
00263       shlon(0) =    0.0
00264       shlon(1) = -360.0
00265       shlon(2) =  360.0
00266 !
00267       n_corners = nca (interpolation_mode)
00268 !
00269       length (1:ndim_3d) = grid_shape(2,1:ndim_3d) - &
00270                            grid_shape(1,1:ndim_3d) + 1
00271 !
00272       grid_id = use_how(3)
00273       gp => Grids(grid_id)
00274 !
00275 ! ndim = ndim_2d: Locations for 2D conservative remapping
00276 ! ndim = ndim_3d: Locations for 3D conservative remapping
00277 !
00278       ndim = interpolation_mode
00279 !
00280 !===> Find overlapping cells
00281 !
00282       nprev = 0
00283 
00284       nn = 0
00285 
00286       allocate(tgt(4), src(4), STAT=ierror)
00287 
00288       if ( ierror > 0 ) then
00289          ierrp (1) = ierror
00290          ierrp (2) = 8
00291          ierror = PRISM_Error_Alloc
00292          call psmile_error ( ierror, 'tgt and src struct', &
00293               ierrp, 2, __FILE__, __LINE__ )
00294          return
00295       endif
00296 
00297       allocate(visited(grid_shape(1,1):grid_shape(2,1), &
00298                        grid_shape(1,2):grid_shape(2,2)), STAT=ierror)
00299 
00300       if ( ierror > 0 ) then
00301          ierrp (1) = ierror
00302          ierrp (2) = length(1) * length(2)
00303          ierror = PRISM_Error_Alloc
00304          call psmile_error ( ierror, 'visited', &
00305               ierrp, 2, __FILE__, __LINE__ )
00306          return
00307       endif
00308 
00309       visited   = 0
00310 
00311       bnd_size = bnd_size_inc
00312 
00313       allocate(search%boundary_cell(bnd_size,5), STAT=ierror)
00314 
00315       if ( ierror > 0 ) then
00316          ierrp (1) = ierror
00317          ierrp (2) = 5*bnd_size
00318          ierror = PRISM_Error_Alloc
00319          call psmile_error ( ierror, 'boundary_cell', &
00320               ierrp, 2, __FILE__, __LINE__ )
00321          return
00322       endif
00323 
00324       search%boundary_cell(:,:) = PSMILe_Undef
00325 
00326       stacksize = stacksize_inc
00327 
00328       allocate(stack(stacksize), STAT=ierror)
00329 
00330       if ( ierror > 0 ) then
00331          ierrp (1) = ierror
00332          ierrp (2) = stacksize
00333          ierror = PRISM_Error_Alloc
00334          call psmile_error ( ierror, 'stack', &
00335               ierrp, 2, __FILE__, __LINE__ )
00336          return
00337       endif
00338 
00339       nbr_cells_tot = nbr_cells_inc
00340 
00341       allocate(tmp_neighcells_3d(3,nbr_cells_tot), STAT=ierror)
00342 
00343       if ( ierror > 0 ) then
00344          ierrp (1) = ierror
00345          ierrp (2) = 2 * nbr_cells_tot
00346          ierror = PRISM_Error_Alloc
00347          call psmile_error ( ierror, 'tmp_neighcells_3d', &
00348               ierrp, 2, __FILE__, __LINE__ )
00349          return
00350       endif
00351 
00352       do ipart = 1, search%npart
00353 
00354          ntgtlocs(1) = count(msklocs(1,ipart)%vector)
00355          ntgtlocs(2) = count(msklocs(2,ipart)%vector)
00356 
00357 !  lenij  = Total number of entries in first 2 directions
00358 !  lenijk = Total number of entries in all   3 directions
00359 
00360          lenij  = ntgtlocs (1)
00361          lenijk = ntgtlocs (1) * ntgtlocs (2)
00362 
00363          nend = nprev + lenijk
00364 
00365          if ( lenijk == 0 ) cycle
00366 
00367 #ifdef PRISM_ASSERTION
00368 !
00369 !===> Is the number of neighbors (number of corners) sufficient
00370 !
00371          if (num_neigh < n_corners) then
00372             print 9970, trim(ch_id), num_neigh, n_corners
00373             call psmile_assert (__FILE__, __LINE__, &
00374                  'Number of neighbors is insufficient')
00375          endif
00376 !
00377 !===> 
00378 !
00379          if (ncpl < nprev + ntgtlocs (1) * ntgtlocs (2) ) then
00380             print *, trim(ch_id), " ncpl, nprev, ntgtlocs ", ncpl, nprev, ntgtlocs
00381             call psmile_assert (__FILE__, __LINE__, &
00382                  'ncpl < nprev + PRODUCT (ntgtlocs) ')
00383          endif
00384 !
00385 !===> Are the locations within the correct shape ?
00386 !     Note: the scrloc's are source locations in the 
00387 !           method grid which has an virtual cell 0.
00388 !
00389 !cdir vector
00390          if ( use_how(2) == PRISM_Gaussreduced_regvrt ) then
00391             do i = 1, npoints(1,ipart)
00392                if ( abs(srclocs(1,ipart)%vector(2*i-1)) < grid_shape(1,1) .or. &
00393                     abs(srclocs(1,ipart)%vector(2*i-1)) > grid_shape(2,1) .or. &
00394                     abs(srclocs(1,ipart)%vector(2*i  )) < grid_shape(1,2) .or. &
00395                     abs(srclocs(1,ipart)%vector(2*i  )) > grid_shape(2,2)) exit
00396             end do
00397          else
00398             do i = 1, npoints(1,ipart)
00399                if ( srclocs(1,ipart)%vector(2*i-1) < grid_shape(1,1) .or. &
00400                     srclocs(1,ipart)%vector(2*i-1) > grid_shape(2,1) .or. &
00401                     srclocs(1,ipart)%vector(2*i  ) < grid_shape(1,2) .or. &
00402                     srclocs(1,ipart)%vector(2*i  ) > grid_shape(2,2)) exit
00403             end do
00404          endif
00405 
00406          if (i <= npoints(1,ipart)) then
00407             print *, "Incorrect index in srclocs, i", &
00408                  i, srclocs(1,ipart)%vector(2*i-1), srclocs(1,ipart)%vector(2*i)
00409             call psmile_assert (__FILE__, __LINE__, 'wrong index')
00410          endif
00411 
00412 !cdir vector
00413          do i = 1, npoints(2,ipart)
00414             if ( srclocs(2,ipart)%vector(i) < grid_shape(1,3)-1 .or. &
00415                  srclocs(2,ipart)%vector(i) > grid_shape(2,3)) exit
00416          end do
00417 
00418          if (i <= npoints(2,ipart)) then
00419             print *, "Incorrect index in srclocs(2), i", i, srclocs(2,ipart)%vector(i), &
00420                  grid_shape(1,3)-1, grid_shape(2,3)
00421             call psmile_assert (__FILE__, __LINE__, 'wrong index')
00422          endif
00423 #endif
00424          m    = 0
00425          mm   = 0
00426          mold = 0
00427 
00428          do n = 1, npoints(1,ipart)
00429 
00430             ! msklocs contains source cells that shall not be considered
00431             ! as starting cells for further local search. In particular
00432             ! this is required for reglonlatvrt source cell grids.
00433 
00434             if ( .not. msklocs(1,ipart)%vector(n) ) cycle
00435 
00436             m  = m + 1
00437 
00438             do j1 = 0, tgt_corners-1
00439                j2 = mod(j1+1,tgt_corners)
00440                tgt(j1+1)%p1%x = tgt_cell(1)%vector(j1*ncpl+nprev+m)
00441                tgt(j1+1)%p1%y = tgt_cell(2)%vector(j1*ncpl+nprev+m)
00442                tgt(j1+1)%p2%x = tgt_cell(1)%vector(j2*ncpl+nprev+m)
00443                tgt(j1+1)%p2%y = tgt_cell(2)%vector(j2*ncpl+nprev+m)
00444             enddo
00445 
00446             tgt_min =  999.0d0
00447             tgt_max = -999.0d0
00448 
00449             do j1 = 1, tgt_corners
00450               tgt_min = min(tgt_min,tgt(j1)%p1%x)
00451               tgt_max = min(tgt_max,tgt(j1)%p1%x)
00452             enddo
00453             !
00454             ! convert target longitudes to [0,360] interval if necessary
00455             !
00456             if ( tgt_max - tgt_min > 180.0d0 ) then
00457 
00458                do j1 = 1, tgt_corners
00459 
00460                   if ( tgt(j1)%p1%x > 360.0d0 ) then
00461                       tgt(j1)%p1%x = tgt(j1)%p1%x - 360.0d0
00462                       tgt_cell(1)%vector((j1-1)*ncpl+nprev+m) = &
00463                           tgt_cell(1)%vector((j1-1)*ncpl+nprev+m) - 360.0d0
00464                   endif
00465 
00466                   if ( tgt(j1)%p1%x < zero    ) then
00467                       tgt(j1)%p1%x = tgt(j1)%p1%x + 360.0d0
00468                       tgt_cell(1)%vector((j1-1)*ncpl+nprev+m) = &
00469                           tgt_cell(1)%vector((j1-1)*ncpl+nprev+m) + 360.0d0
00470                   endif
00471 
00472                   if ( tgt(j1+1)%p2%x > 360.0d0 ) then
00473                       tgt(j1+1)%p2%x = tgt(j1+1)%p2%x - 360.0d0
00474                       j2 = mod(j1,tgt_corners)
00475                       tgt_cell(1)%vector(j2*ncpl+nprev+m) = &
00476                               tgt_cell(1)%vector(j2*ncpl+nprev+m) - 360.0d0
00477                   endif
00478 
00479                   if ( tgt(j1+1)%p2%x < zero    ) then
00480                       tgt(j1+1)%p2%x = tgt(j1+1)%p2%x + 360.0d0
00481                       tgt_cell(1)%vector(j2*ncpl+nprev+m) = &
00482                               tgt_cell(1)%vector(j2*ncpl+nprev+m) + 360.0d0
00483                   endif
00484 
00485                enddo
00486 
00487             endif
00488 
00489             !
00490             ! For each new target cell the stack is reused and
00491             ! the first entry in the stack is initialised.
00492             !
00493             index = 1
00494 
00495             stack(index)%i   = srclocs(1,ipart)%vector(2*n-1)
00496             stack(index)%j   = srclocs(1,ipart)%vector(2*n  )
00497             stack(index)%dir = 0
00498             !
00499             ! ... maybe we should first test whether our initial guess is really good
00500             !
00501             visited(stack(index)%i,stack(index)%j) = nprev+m
00502             !
00503             ! Store results of first cell
00504             !
00505             ! ... allocate new memory if necessary
00506             !
00507             if ( nn + 1 > nbr_cells_tot ) then
00508 
00509                nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00510 
00511                allocate(new_neighcells_3d(2,nbr_cells_tot), STAT=ierror)
00512 
00513                if ( ierror > 0 ) then
00514                   ierrp (1) = ierror
00515                   ierrp (2) = 2 * nbr_cells_tot
00516                   ierror = PRISM_Error_Alloc
00517                   call psmile_error ( ierror, 'new_neighcells_3d', &
00518                        ierrp, 2, __FILE__, __LINE__ )
00519                   return
00520                endif
00521 
00522                new_neighcells_3d(1:2,1:nn) = tmp_neighcells_3d(1:2,1:nn)
00523 
00524                deallocate(tmp_neighcells_3d, STAT=ierror)
00525 
00526                if ( ierror > 0 ) then
00527                   ierrp (1) = ierror
00528                   ierrp (2) = 3 * ( nbr_cells_tot - nbr_cells_inc )
00529                   ierror = PRISM_Error_Dealloc
00530                   call psmile_error ( ierror, 'tmp_neighcells_3d', &
00531                        ierrp, 2, __FILE__, __LINE__ )
00532                   return
00533                endif
00534 
00535                tmp_neighcells_3d => new_neighcells_3d
00536 
00537             endif
00538 
00539             nn = nn + 1
00540 
00541             tmp_neighcells_3d(1,nn) = stack(index)%i
00542             tmp_neighcells_3d(2,nn) = stack(index)%j
00543 
00544             nbr_cells(nprev+m) = 1
00545 
00546 #ifdef META_POST
00547             if ( nprev+m == GRIDPOINT ) then
00548 
00549                ! 595 points width for an A4 page, scaled by 180
00550                ! a shift my have to be applied in order to get
00551                ! the boxes on display
00552                
00553                shift_x = -tgt(1)%p1%x+180.0/6.0
00554                shift_y = -tgt(1)%p1%y+180.0/6.0
00555 
00556                open ( unit=99, file='cell_search.txt', form='formatted')
00557                write(99,'(a)') 'beginfig(-1)'
00558                write(99,'(a)') 'numeric u;'
00559                write(99,'(a3,f8.4,a1)') 'u:=', 595.0/180, ';'
00560                write(99,'(a)') 'pickup pencircle scaled .1pt;'
00561 
00562                ii = stack(index)%i
00563                jj = stack(index)%j
00564 
00565                if ( gp%grid_type == PRISM_Reglonlatvrt ) then
00566                   do j1 = 1, n_corners
00567                      j2 = mod(j1,tgt_corners)+1
00568                      src(j1)%p1%x = corner_x(ii, 1,edge(1,j1))
00569                      src(j1)%p1%y = corner_y( 1,jj,edge(2,j1))
00570                      src(j1)%p2%x = corner_x(ii, 1,edge(1,j2))
00571                      src(j1)%p2%y = corner_y( 1,jj,edge(2,j2))
00572                   enddo
00573                else
00574                   do j1 = 1, n_corners
00575                      j2 = mod(j1,tgt_corners)+1
00576                      src(j1)%p1%x = corner_x(ii,jj,j1)
00577                      src(j1)%p1%y = corner_y(ii,jj,j1)
00578                      src(j1)%p2%x = corner_x(ii,jj,j2)
00579                      src(j1)%p2%y = corner_y(ii,jj,j2)
00580                   enddo
00581                endif
00582 
00583                ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00584 
00585                do j1 = 1, n_corners
00586                   print *, ' source cell point o', ii, jj, n, src(j1)%p1%x, src(j1)%p1%y
00587                   ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00588                enddo
00589                print *, ' --------- '
00590 
00591             endif
00592 #endif
00593             !
00594             ! Store further results with Andreas algorithm
00595             !
00596             unfolded_check = 0
00597 
00598             do while ( index >= 0 )
00599 
00600                if ( index == 0 .and. unfolded_check > 0 ) then
00601 
00602                   call psmile_unfolded_check_dble( m,  & 
00603                                       unfolded_check,  &
00604                                       gp%grid_type,    &
00605                                       visited,         & 
00606                                       tgt_corners,     & 
00607                                       tgt,             & 
00608                                       grid_shape,      & 
00609                                       corner_shape_3d, & 
00610                                       n_corners,       & 
00611                                       corner_x,        & 
00612                                       corner_y,        & 
00613                                       stacksize,       & 
00614                                       stack,           & 
00615                                       index,           & 
00616                                       ierror )
00617 
00618                   unfolded_check = 0
00619 
00620                   if ( index > 0 ) then
00621                      !
00622                      ! Store results of first cell
00623                      !
00624                      ! ... allocate new memory if necessary
00625                      !
00626                      if ( nn + 1 > nbr_cells_tot ) then
00627 
00628                         nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00629 
00630                         allocate(new_neighcells_3d(2,nbr_cells_tot), STAT=ierror)
00631 
00632                         if ( ierror > 0 ) then
00633                            ierrp (1) = ierror
00634                            ierrp (2) = 2 * nbr_cells_tot
00635                            ierror = PRISM_Error_Alloc
00636                            call psmile_error ( ierror, 'new_neighcells_3d', &
00637                                 ierrp, 2, __FILE__, __LINE__ )
00638                            return
00639                         endif
00640 
00641                         new_neighcells_3d(1:2,1:nn) = tmp_neighcells_3d(1:2,1:nn)
00642 
00643                         deallocate(tmp_neighcells_3d, STAT=ierror)
00644 
00645                         if ( ierror > 0 ) then
00646                            ierrp (1) = ierror
00647                            ierrp (2) = 2 * ( nbr_cells_tot - nbr_cells_inc )
00648                            ierror = PRISM_Error_Dealloc
00649                            call psmile_error ( ierror, 'tmp_neighcells_3d', &
00650                                 ierrp, 2, __FILE__, __LINE__ )
00651                            return
00652                         endif
00653 
00654                         tmp_neighcells_3d => new_neighcells_3d
00655 
00656                      endif
00657 
00658                      nn = nn + 1
00659 
00660                      tmp_neighcells_3d(1,nn) = stack(index)%i
00661                      tmp_neighcells_3d(2,nn) = stack(index)%j
00662 
00663                      nbr_cells(nprev+m) = nbr_cells(nprev+m) + 1
00664 
00665                      cycle
00666 
00667                   endif ! index > 0
00668 
00669                endif ! index == 0 .and. unfolded_check > 0
00670 
00671                if ( index == 0 ) exit ! while loop
00672 
00673                !
00674                ! Change direction
00675                !
00676                stack(index)%dir = stack(index)%dir + 1
00677 
00678                if ( stack(index)%dir > 4 ) then
00679                   !
00680                   ! All 4 directions are already tested
00681                   ! for this cell so we have to go back
00682                   !
00683                   index = index - 1
00684                   cycle
00685                endif
00686                !
00687                ! Go to the next cell
00688                !
00689                ii = stack(index)%i + ijdst(1,stack(index)%dir)
00690                jj = stack(index)%j + ijdst(2,stack(index)%dir)
00691                !
00692                ! Apply cyclic boundary conditions if required
00693                !
00694                if ( cyclic(1) ) ii = mod(ii+length(1)-grid_shape(1,1),length(1))+grid_shape(1,1)
00695                if ( cyclic(2) ) jj = mod(jj+length(2)-grid_shape(1,2),length(2))+grid_shape(1,2)
00696                !
00697                ! Check that indices are in the allowed range
00698                !
00699                if ( ii < grid_shape(1,1) .or. ii > grid_shape(2,1) .or. &
00700                     jj < grid_shape(1,2) .or. jj > grid_shape(2,2) ) then
00701 
00702                   if ( gp%grid_type == PRISM_Irrlonlat_regvrt ) then
00703                      if ( jj < grid_shape(1,2) ) unfolded_check = 1
00704                      if ( jj > grid_shape(2,2) ) unfolded_check = 2
00705                   endif
00706 
00707 
00708                   if ( m > mold ) then
00709 
00710                      mold = m
00711 
00712                      if ( mm + 1 > bnd_size ) then
00713 
00714                         bnd_size = bnd_size + bnd_size_inc
00715 
00716                         allocate(boundary_cell(bnd_size,5), STAT=ierror)
00717 
00718                         if ( ierror > 0 ) then
00719                            ierrp (1) = ierror
00720                            ierrp (2) = 5 * bnd_size
00721                            ierror = PRISM_Error_Alloc
00722                            call psmile_error ( ierror, 'boundary_cell', &
00723                                 ierrp, 2, __FILE__, __LINE__ )
00724                            return
00725                         endif
00726 
00727                         boundary_cell(1:mm,1:5) = search%boundary_cell(1:mm,1:5)
00728 
00729                         boundary_cell(mm+1:bnd_size,1:5) = PSMILe_Undef
00730 
00731                         deallocate(search%boundary_cell, STAT=ierror)
00732 
00733                         if ( ierror > 0 ) then
00734                            ierrp (1) = ierror
00735                            ierrp (2) = 5 * ( bnd_size - bnd_size_inc )
00736                            ierror = PRISM_Error_Dealloc
00737                            call psmile_error ( ierror, 'boundary_cell', &
00738                                 ierrp, 2, __FILE__, __LINE__ )
00739                            return
00740                         endif
00741 
00742                         search%boundary_cell => boundary_cell
00743 
00744                      endif
00745 
00746                      mm = mm + 1
00747 
00748                      search%boundary_cell(mm,1) = m
00749                      search%boundary_cell(mm,5) = psmile_rank
00750 
00751                   endif
00752 
00753                   cycle
00754 
00755                endif
00756                !
00757                ! When cell is already checked for this target ...
00758                !
00759                if ( visited(ii,jj) == nprev+m ) cycle
00760                !
00761                ! Test current cell
00762                !
00763                if ( gp%grid_type == PRISM_Reglonlatvrt ) then
00764                   do j1 = 1, n_corners
00765                      j2 = mod(j1,tgt_corners)+1
00766                      src(j1)%p1%x = corner_x(ii, 1,edge(1,j1))
00767                      src(j1)%p1%y = corner_y( 1,jj,edge(2,j1))
00768                      src(j1)%p2%x = corner_x(ii, 1,edge(1,j2))
00769                      src(j1)%p2%y = corner_y( 1,jj,edge(2,j2))
00770                   enddo
00771                else
00772                   do j1 = 1, n_corners
00773                      j2 = mod(j1,tgt_corners)+1
00774                      src(j1)%p1%x = corner_x(ii,jj,j1)
00775                      src(j1)%p1%y = corner_y(ii,jj,j1)
00776                      src(j1)%p2%x = corner_x(ii,jj,j2)
00777                      src(j1)%p2%y = corner_y(ii,jj,j2)
00778                   enddo
00779                endif
00780 
00781                call psmile_overlap_dble ( n_corners, tgt_corners, src, tgt, overlap )
00782                !
00783                ! Mark cell as checked for this cycle
00784                !
00785                visited(ii,jj) = nprev+m
00786 
00787 #ifdef META_POST
00788                if ( nprev+m == GRIDPOINT ) then
00789 
00790                   if ( overlap ) then
00791                      ierror = mp_fill_ov ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00792                      do j1 = 1,  n_corners
00793                         print *, ' source cell point o', ii, jj, n, src(j1)%p1%x, src(j1)%p1%y
00794                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00795                      enddo
00796                   else
00797                      ierror = mp_fill_vi ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1 )
00798                      do j1 = 1,  n_corners
00799                         print *, ' source cell point v', ii, jj, n, src(j1)%p1%x, src(j1)%p1%y
00800                         ierror = mp_draw_red ( src(j1)%p1, src(j1)%p2 )
00801                      enddo
00802                   endif
00803                   print *, ' --------- '
00804 
00805                endif
00806 #endif
00807                if ( .not. overlap ) then
00808 
00809                   cycle
00810 
00811                else
00812                   !
00813                   ! We found a new source cell
00814                   !
00815                   ! ... allocate new memory if necessary
00816                   !
00817                   if ( nn + 1 > nbr_cells_tot ) then
00818 
00819                      nbr_cells_tot = nbr_cells_tot + nbr_cells_inc
00820 
00821                      allocate(new_neighcells_3d(2,nbr_cells_tot), STAT=ierror)
00822 
00823                      if ( ierror > 0 ) then
00824                         ierrp (1) = ierror
00825                         ierrp (2) = 2 * nbr_cells_tot
00826                         ierror = PRISM_Error_Alloc
00827                         call psmile_error ( ierror, 'new_neighcells_3d', &
00828                              ierrp, 2, __FILE__, __LINE__ )
00829                         return
00830                      endif
00831 
00832                      new_neighcells_3d(1:2,1:nn) = tmp_neighcells_3d(1:2,1:nn)
00833 
00834                      deallocate(tmp_neighcells_3d, STAT=ierror)
00835 
00836                      if ( ierror > 0 ) then
00837                         ierrp (1) = ierror
00838                         ierrp (2) = 2 * ( nbr_cells_tot - nbr_cells_inc )
00839                         ierror = PRISM_Error_Dealloc
00840                         call psmile_error ( ierror, 'tmp_neighcells_3d', &
00841                              ierrp, 2, __FILE__, __LINE__ )
00842                         return
00843                      endif
00844 
00845                      tmp_neighcells_3d => new_neighcells_3d
00846 
00847                   endif
00848                   !
00849                   ! Store results
00850                   !
00851                   nn = nn + 1
00852 
00853                   tmp_neighcells_3d(1,nn) = ii
00854                   tmp_neighcells_3d(2,nn) = jj
00855 
00856                   nbr_cells(nprev+m) = nbr_cells(nprev+m) + 1
00857 
00858                   !
00859                   ! Prepare for the next source cell to test
00860                   !
00861                   ! ... allocate new memory for the stack if necessary
00862                   !
00863                   if ( index + 1 > stacksize ) then
00864                      stacksize =  stacksize + stacksize_inc
00865 
00866                      allocate(new_stack(stacksize), STAT=ierror)
00867 
00868                      if ( ierror > 0 ) then
00869                         ierrp (1) = ierror
00870                         ierrp (2) = stacksize
00871                         ierror = PRISM_Error_Alloc
00872                         call psmile_error ( ierror, 'stack', &
00873                              ierrp, 2, __FILE__, __LINE__ )
00874                         return
00875                      endif
00876 
00877                      new_stack(1:index) = stack(1:index)
00878 
00879                      deallocate(stack, STAT=ierror)
00880 
00881                      if ( ierror > 0 ) then
00882                         ierrp (1) = ierror
00883                         ierrp (2) = stacksize - stacksize_inc
00884                         ierror = PRISM_Error_Dealloc
00885                         call psmile_error ( ierror, 'stack', &
00886                              ierrp, 2, __FILE__, __LINE__ )
00887                         return
00888                      endif
00889 
00890                      stack => new_stack
00891 
00892                   endif
00893                   !
00894                   ! ... initialise new stack entry
00895                   !
00896                   index = index + 1
00897                   stack(index)%i   = ii
00898                   stack(index)%j   = jj
00899                   stack(index)%dir = 0
00900 
00901                endif
00902 
00903             enddo ! index
00904 #ifdef META_POST
00905             if ( nprev+m == GRIDPOINT ) then
00906                do j1 = 1, tgt_corners
00907                   ierror = mp_draw_black ( tgt(j1)%p1, tgt(j1)%p2 )
00908                   print *, ' target cell point ', nn, n, tgt(j1)%p1%x, tgt(j1)%p1%y
00909                enddo
00910 
00911                write(99,'(a)') 'endfig'
00912                write(99,'(a)') 'end'
00913                close(unit=99)
00914 
00915             endif
00916             ! for later graphical analysis print number of cells found for each target
00917             print *, ' nbr_cells(m) ', m, nbr_cells(m)
00918 #endif
00919          enddo ! npoints
00920          !
00921          ! Apply the results for the first level to all remaining levels
00922          !
00923          do k = 2, ntgtlocs (2)
00924             nbr_cells(nprev+(k-1)*ntgtlocs(1)+1:nprev+k*ntgtlocs(1)) = nbr_cells(nprev+1:nprev+ntgtlocs(1))
00925          end do
00926 
00927          nprev = nend
00928 
00929       enddo ! ipart
00930 
00931       nbr_cells_tot = sum(nbr_cells)
00932       !
00933       ! Allocate memory for final results
00934       !
00935       allocate(neighcells_3d(3,nbr_cells_tot,1), STAT=ierror)
00936 
00937       if ( ierror > 0 ) then
00938          ierrp (1) = ierror
00939          ierrp (2) = 3 * nbr_cells_tot
00940          ierror = PRISM_Error_Alloc
00941          call psmile_error ( ierror, 'neighcells_3d', &
00942               ierrp, 2, __FILE__, __LINE__ )
00943          return
00944       endif
00945       !
00946       ! Copy temporary results to final destination
00947       !
00948       offset1 = 0
00949       offset2 = 0
00950       nprev   = 0
00951 
00952       do ipart = 1, search%npart
00953 
00954          ntgtlocs(1) = count(msklocs(1,ipart)%vector)
00955          ntgtlocs(2) = count(msklocs(2,ipart)%vector)
00956 
00957          if ( ntgtlocs(1) * ntgtlocs(2) == 0 ) cycle
00958 
00959          !  lenij  = Total number of entries in first 2 directions
00960          !  lenijk = Total number of entries in all   3 directions
00961 
00962          lenij  = sum(nbr_cells(offset2+1:offset2+ntgtlocs(1)))
00963          lenijk = lenij * ntgtlocs (2)
00964 
00965          nend = nprev + lenij
00966 
00967          neighcells_3d (1,offset1+1:offset1+lenij,1) = tmp_neighcells_3d(1,nprev+1:nend)
00968          neighcells_3d (2,offset1+1:offset1+lenij,1) = tmp_neighcells_3d(2,nprev+1:nend)
00969 
00970          ! Note: At this point offset contains the previous 3d part plus the first k-level
00971          !       of the current part.
00972 
00973          do k = 2, ntgtlocs (2)
00974             neighcells_3d (1,offset1+(k-1)*lenij+1: offset1+k*lenij,1) = & 
00975                  neighcells_3d (1,offset1+1:offset1+1+lenij,1)
00976             neighcells_3d (2,offset1+(k-1)*lenij+1: offset1+k*lenij,1) = &
00977                  neighcells_3d (2,offset1+1:offset1+1+lenij,1)
00978          end do
00979 
00980          do k = 1, ntgtlocs(2)
00981             neighcells_3d (3,offset1+(k-1)*lenij+1: offset1+k*lenij,1) = &
00982                  srclocs(2,ipart)%vector(k) - grid_shape(1,3) + 1
00983          enddo
00984 
00985          offset1 = offset1 + lenijk
00986          offset2 = offset2 + ntgtlocs(1)*ntgtlocs(2)
00987 
00988          nprev = nend
00989 
00990       enddo ! ipart
00991 !
00992 !===> Control special cases (2d hyperplanes)
00993 !
00994 
00995 !
00996 !===> Deallocate local memory
00997 !
00998       deallocate(tmp_neighcells_3d, STAT=ierror)
00999 
01000       if ( ierror > 0 ) then
01001          ierrp (1) = ierror
01002          ierrp (2) = 2 * nbr_cells_tot
01003          ierror = PRISM_Error_Dealloc
01004          call psmile_error ( ierror, 'tmp_neighcells_3d', &
01005               ierrp, 2, __FILE__, __LINE__ )
01006          return
01007       endif
01008 
01009       deallocate(stack, STAT=ierror)
01010 
01011       if ( ierror > 0 ) then
01012          ierrp (1) = ierror
01013          ierrp (2) = stacksize
01014          ierror = PRISM_Error_Dealloc
01015          call psmile_error ( ierror, 'stack', &
01016               ierrp, 2, __FILE__, __LINE__ )
01017          return
01018       endif
01019 !
01020 !===> All done
01021 !
01022 #ifdef VERBOSE
01023       print 9980, trim(ch_id)
01024 
01025       call psmile_flushstd
01026 #endif /* VERBOSE */
01027 !
01028 !  Formats:
01029 !
01030 9990 format (1x, a, ': psmile_neigh_cells_irreg2_dble')
01031 9980 format (1x, a, ': psmile_neigh_cells_irreg2_dble: eof')
01032 
01033 #ifdef PRISM_ASSERTION
01034 9970 format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
01035 #endif
01036 
01037       end subroutine psmile_neigh_cells_irreg2_dble