psmile_neigh_cells_3d_reg_real.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2007-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_Neigh_cells_3d_reg_real
00008 !
00009 ! !INTERFACE:
00010       subroutine psmile_neigh_cells_3d_reg_real (    &
00011                grid_valid_shape, interpolation_mode, &
00012                cyclic, grid_id, search, corners,     &
00013                npoints, srclocs, ncpl, nbr_cells,    &
00014                ierror )
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019 !
00020       use PSMILe, dummy_interface => PSMILe_Neigh_cells_3d_reg_real
00021 
00022       implicit none
00023 !
00024 ! !INPUT PARAMETERS:
00025 !
00026       Integer, Intent (In)               :: grid_valid_shape (2, ndim_3d)
00027 
00028 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00029 !
00030       Integer,            Intent (In)    :: interpolation_mode
00031 !
00032 !     Specifies the interpolation mode with
00033 !        interpolation_mode = 3 : 3D conservative remapping
00034 !                                 (determines 8 corners)
00035 !        interpolation_mode = 2 : 2D conservative remapping
00036 !                                 (determines 4 corners)
00037 !        interpolation_mode = 1 : 1D conservative remapping
00038 !                                 (determines 2 corners)
00039 !
00040       Logical,            Intent (In)    :: cyclic (ndim_3d)
00041 !    
00042 !     Does the block have cyclic coordinates in the corresponding direction ?
00043 !
00044       Integer,            Intent (In)    :: grid_id
00045 
00046 !     Grid Id of the (local) donor grid
00047  
00048       Type (Enddef_search), Intent (In)  :: search
00049 
00050 !     Info's on coordinates to be searched
00051 !
00052       Type (real_vector), Intent (In)    :: corners (ndim_3d, search%npart)
00053 
00054 !     Coordinates to be searched
00055 
00056       Integer, Intent (In)               :: npoints (ndim_3d, search%npart)
00057 
00058 !     Total number of locations to be transferred
00059 
00060       Type (integer_vector), Intent (In) :: srclocs (ndim_3d, search%npart)
00061 
00062 !     Indices of the grid cell in which the point was found.
00063 !
00064       Integer, Intent (In)               :: ncpl
00065 !
00066 !     Total number of locations to be transferred
00067 !
00068 ! !OUTPUT PARAMETERS:
00069 !
00070       Integer, Intent (Out)              :: nbr_cells(ncpl)
00071 !
00072 !     Number of cells found per target cell
00073 !
00074       Integer, Intent (Out)              :: ierror
00075 
00076 !     Returns the error code of PSMILE_Neigh_cells_3d_reg_real;
00077 !             ierror = 0 : No error
00078 !             ierror > 0 : Severe error
00079 !
00080 ! !DEFINED PARAMETERS:
00081 !
00082 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00083 !               = 2 ** ndim_3d
00084 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00085 !               = 2 ** ndim_2d
00086 !
00087       Integer, Parameter              :: n_corners_3d = 2 * 2 * 2
00088       Integer, Parameter              :: n_corners_2d = 2 * 2
00089       Integer, Parameter              :: n_corners_1d = 2
00090 !
00091 ! !LOCAL VARIABLES
00092 !
00093 !     ... for corner locations searched
00094 !
00095       Real                            :: ctr_corner, global(ndim_3d)
00096       Integer                         :: is
00097 
00098       Type (Corner_Block), Pointer    :: cp
00099 
00100       Type (Dble_vector)              :: chmin(ndim_3d), chmax(ndim_3d)
00101 
00102       Type (Integer_vector)           :: ic (ndim_3d, search%npart)
00103 
00104       Integer                         :: nbr_cells_tot
00105       Integer                         :: ipart, jpart, jstart
00106       Integer                         :: i, j, k, n, nbeg, nend, nprev
00107       Integer                         :: lenij, lenk, lenijk, n_corners
00108       Integer                         :: ii, jj, kk, offset, icell
00109 !
00110 !     ... For locations controlled
00111 !
00112       Integer                         :: length (ndim_3d)
00113       Integer                         :: lenc   (ndim_3d)
00114 !
00115       Integer                         :: nlocs (ndim_3d)
00116       Integer                         :: nca  (ndim_3d)
00117       Integer                         :: ndim ! nbr of dims for conservative remapping
00118       Logical                         :: overlap
00119 !
00120 !     ... for error handling
00121 !
00122       Integer, Parameter              :: nerrp = 1
00123       Integer                         :: ierrp (nerrp)
00124 !
00125 !
00126 ! !DESCRIPTION:
00127 !
00128 ! Subroutine "PSMILe_Neigh_cells_3d_reg_real" searches the cells on the donor
00129 ! grid that overlap with a given target cell.
00130 !
00131 ! Subroutine "PSMILe_Neigh_cells_3d_reg_real" is designed for (target) grids
00132 ! of type "PRISM_reg_reallonlatvrt".
00133 !
00134 ! !TO DO:
00135 !
00136 !===> Check what we have to do for linear interpolation (interpolation_mode)
00137 !      in the vertical. Is this already covered by a subsequent call to
00138 !      neigh_trili_3d?
00139 !
00140 !===> Cyclic boundary conditions
00141 !
00142 !
00143 ! !REVISION HISTORY:
00144 !
00145 !   Date      Programmer   Description
00146 ! ----------  ----------   -----------
00147 ! 06.06.19    R. Redler    created
00148 !
00149 !EOP
00150 !----------------------------------------------------------------------
00151 !
00152 !  $Id: psmile_neigh_cells_3d_reg_real.F90 2082 2009-10-21 13:31:19Z hanke $
00153 !  $Author: hanke $
00154 !
00155    Character(len=len_cvs_string), save :: mycvs = 
00156        '$Id: psmile_neigh_cells_3d_reg_real.F90 2082 2009-10-21 13:31:19Z hanke $'
00157 !
00158 !----------------------------------------------------------------------
00159 !
00160       data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00161       data global /360.0d0, 360.0d0, 1.0d35/
00162 !
00163 !----------------------------------------------------------------------
00164 !
00165 #ifdef VERBOSE
00166       print 9990, trim(ch_id)
00167 
00168       call psmile_flushstd
00169 #endif /* VERBOSE */
00170 !
00171 !===> Initialization
00172 !
00173       ierror = 0
00174 !
00175       cp    => Grids(grid_id)%corner_pointer
00176 !
00177       nprev = 0
00178 
00179       n_corners = nca (interpolation_mode)
00180 !
00181       length (1:ndim_3d) = grid_valid_shape(2,1:ndim_3d) - &
00182                            grid_valid_shape(1,1:ndim_3d) + 1
00183 
00184       lenc (1:ndim_3d)   = cp%corner_shape(2,1:ndim_3d) - &
00185                            cp%corner_shape(1,1:ndim_3d) + 1
00186 !
00187 !===> Control interpolation mode
00188 !
00189       if (interpolation_mode < ndim_2d .or. interpolation_mode > ndim_3d) then
00190          ierrp (1) = interpolation_mode
00191          ierror = PRISM_Error_Internal
00192 
00193          call psmile_error ( ierror, 'unsupported interpolation mode', &
00194                              ierrp, 1, __FILE__, __LINE__ )
00195          return
00196       endif
00197 !
00198 ! ndim = ndim_2d: Locations for 2D conservative remapping
00199 ! ndim = ndim_3d: Locations for 3D conservative remapping
00200 !
00201       ndim = interpolation_mode
00202 
00203 !
00204 !===> Determine extent of each cell 
00205 !
00206       do i = 1, ndim
00207          Allocate( chmin(i)%vector(length(i)), chmax(i)%vector(length(i)), STAT=ierror )
00208          if (ierror /= 0) then
00209             ierrp (1) = 2 * length(i)
00210             ierror = PRISM_Error_Alloc
00211             call psmile_error ( ierror, 'chmin, chmax', &
00212                  ierrp, 1, __FILE__, __LINE__ )
00213             return
00214          endif
00215       enddo
00216 
00217       do i = 1, ndim
00218          jstart = Grids(grid_id)%grid_shape(1,i) - cp%corner_shape(1,i)
00219          do jpart = 1, length(i), 256
00220 !cdir vector
00221             do j = jpart, min(length(i),jpart+256-1)
00222                chmin(i)%vector(j) = min(cp%corners_real(i)%vector(jstart+j), &
00223                                         cp%corners_real(i)%vector(jstart+lenc(i)+j))
00224                chmax(i)%vector(j) = max(cp%corners_real(i)%vector(jstart+j), &
00225                                         cp%corners_real(i)%vector(jstart+lenc(i)+j))
00226             enddo
00227          enddo
00228       enddo
00229 
00230       do ipart = 1, search%npart
00231          do i = 1, ndim
00232             Allocate( ic(i,ipart)%vector(npoints(i,ipart)-1), STAT=ierror )
00233             if (ierror /= 0) then
00234                ierrp (1) = npoints(i,ipart)
00235                ierror = PRISM_Error_Alloc
00236                call psmile_error ( ierror, 'ic vector', &
00237                     ierrp, 1, __FILE__, __LINE__ )
00238                return
00239             endif
00240          enddo
00241       enddo
00242 
00243       do ipart = 1, search%npart
00244 
00245          nlocs(1) = npoints(1,ipart) - 1
00246          nlocs(2) = npoints(2,ipart) - 1
00247          if ( ndim == ndim_3d ) then
00248             nlocs(3) = npoints(3,ipart) - 1
00249          else
00250             nlocs(3) = npoints(3,ipart)
00251          endif
00252 
00253 !  lenij  = Total number of entries in first 2 directions
00254 !  lenijk = Total number of entries in all   3 directions
00255 
00256          lenij  = nlocs (1) * nlocs (2)
00257          lenijk = lenij * nlocs (3)
00258 
00259          nbeg = nprev + 1
00260          nend = nprev + lenijk
00261 !
00262 #ifdef PRISM_ASSERTION
00263 !
00264 !===> Is the dimension of array "neighbors_3d" sufficient
00265 !
00266          if (ncpl < nprev + lenijk ) then
00267             print *, trim(ch_id), " ncpl, nprev, nlocs ", ncpl, nprev, nlocs
00268             call psmile_assert (__FILE__, __LINE__, &
00269                  'ncpl < nprev + lenijk ')
00270          endif
00271 !
00272 !===> Are the locations within the correct shape ?
00273 !     Note: the scrloc's point to the source cell in the local grid.
00274 !
00275          do j = 1, ndim_3d
00276 !cdir vector
00277             do i = 1, nlocs (j)
00278                if (srclocs(j,ipart)%vector(i) < grid_valid_shape(1,j) .or. &
00279                    srclocs(j,ipart)%vector(i) > grid_valid_shape(2,j))  exit
00280             end do
00281 
00282             if (i <= nlocs(j)) then
00283                print *, trim(ch_id),':','Incorrect index in j, i', &
00284                                          j, i, srclocs(j,ipart)%vector(i)
00285                call psmile_assert (__FILE__, __LINE__, &
00286                     'wrong index')
00287             endif
00288 
00289          enddo
00290 #endif
00291 !
00292 !===> Calculate number of donor cells needed for each target cell
00293 !
00294          do i = 1, ndim
00295             do jpart = 1, nlocs(i), 256
00296 !cdir vector
00297             do j = jpart, min(nlocs(i),jpart+256-1)
00298 
00299                do n = srclocs(i,ipart)%vector(j) - grid_valid_shape(1,i) + 1, length(i)
00300 
00301                  if  ( ( corners(i,ipart)%vector(j+1) - chmin(i)%vector(n) ) * &
00302                        ( chmax(i)%vector(n) - corners(i,ipart)%vector(j) ) < 0 ) &
00303                   exit  ! if no overlap
00304 
00305                enddo
00306 
00307                ! the number of overlapping cells is ...
00308                ic(i,ipart)%vector(j) = n - srclocs(i,ipart)%vector(j) + grid_valid_shape(1,i) - 1
00309                
00310                ! Check whether we need to get more cells on the "left" 
00311                
00312                is = srclocs(i,ipart)%vector(j) - grid_valid_shape(1,i) + 1
00313 
00314                if ( corners(i,ipart)%vector(j) < chmin(i)%vector(is) .and. cyclic(i) ) then
00315 #ifdef DEBUG
00316                   print *,  trim(ch_id),':',' search cell minimum      ', corners(i,ipart)%vector(j)
00317                   print *,  trim(ch_id),':',' is outside of local cell ', chmin(i)%vector(is)
00318 #endif
00319                   n = minval(srclocs(i,ipart)%vector)
00320 
00321                   do ii = length(i), n-grid_valid_shape(1,i)+1, -1
00322                     if ( corners(i,ipart)%vector(j) < chmax(i)%vector(ii)-global(i) ) exit
00323                   enddo
00324 
00325                   ii = length(i) - ii + 1
00326 
00327                   ic(i,ipart)%vector(j) = ic(i,ipart)%vector(j) + ii
00328 #ifdef DEBUG
00329                   print *,  trim(ch_id),':',' ic vector increased from ', &
00330                                                     ic(i,ipart)%vector(j) -ii, &
00331                                             ' to ', ic(i,ipart)%vector(j)
00332 #endif
00333                endif
00334 
00335 #ifdef DEBUG
00336                if ( ic(i,ipart)%vector(j) == 0 ) then
00337                   n = srclocs(i,ipart)%vector(j) - grid_valid_shape(1,i) + 1
00338                   print *,  trim(ch_id),':',' Problem at     ', chmin(i)%vector(n), &
00339                                                                 chmax(i)%vector(n)
00340                   print *,  trim(ch_id),':',' for element    ', corners(i,ipart)%vector(j), &
00341                                                                 corners(i,ipart)%vector(j+1)
00342                   print *,  trim(ch_id),':',' for with index ', i, ipart, j
00343                endif
00344 #endif
00345 
00346             enddo
00347             enddo
00348          enddo
00349 
00350          if ( ndim == ndim_2d ) then
00351 
00352             do j = 1, nlocs (2)
00353                do i = 1, nlocs (1)
00354                   nbr_cells(nprev+(j-1)*nlocs(1)+i) = ic(1,ipart)%vector(i) &
00355                                                     * ic(2,ipart)%vector(j)
00356                enddo
00357             enddo
00358 
00359             do k = 2, nlocs (3)
00360                nbr_cells (nprev+(k-1)*lenij+1:nprev+k*lenij) = nbr_cells(nprev+1:nprev+lenij)
00361             end do
00362 
00363          else
00364 
00365             do k = 1, nlocs (3)
00366                do j = 1, nlocs (2)
00367                   do i = 1, nlocs (1)
00368                      nbr_cells(nprev+(k-1)*lenij+(j-1)*nlocs(1)+i) = &
00369                                                          ic(1,ipart)%vector(i) &
00370                                                        * ic(2,ipart)%vector(j) &
00371                                                        * ic(3,ipart)%vector(k)
00372                   enddo
00373                enddo
00374             enddo
00375 
00376          endif
00377 !
00378 !===> enddo ipart
00379 !
00380          nprev = nend
00381 
00382       enddo
00383 !
00384 !===> Allocate memory to store corner indices
00385 !
00386 !     Reuse of lenij!
00387 !
00388       if ( ndim == ndim_2d ) then
00389          lenij  = 0
00390          do ipart = 1, search%npart
00391             do j = 1,  npoints(2,ipart)-1
00392                do i = 1,  npoints(1,ipart)-1
00393                   lenij = lenij + ic(1,ipart)%vector(i)*ic(2,ipart)%vector(j)
00394                enddo
00395             enddo
00396          enddo
00397          lenk = npoints(3,1)
00398 #ifdef PRISM_ASSERTION
00399          if ( sum(npoints(3,:))/search%npart /= lenk ) then
00400             print *, trim(ch_id), lenk, npoints(3,:)
00401             call psmile_assert (__FILE__, __LINE__, &
00402                  'npoints(3) should be the same for all parts')
00403          endif
00404 #endif
00405          nbr_cells_tot = lenij * lenk
00406       else
00407 
00408          nbr_cells_tot = sum(nbr_cells)
00409 
00410       endif
00411 
00412       Allocate ( neighcells_3d(ndim_3d,nbr_cells_tot,n_corners), STAT=ierror)
00413 
00414       if (ierror /= 0) then
00415          ierrp (1) = 3 * nbr_cells_tot * n_corners
00416          ierror = PRISM_Error_Alloc
00417          call psmile_error ( ierror, 'neighcells_3d', &
00418               ierrp, 1, __FILE__, __LINE__ )
00419          return
00420       endif
00421 
00422       neighcells_3d = PSMILe_undef
00423 
00424       if ( ndim == ndim_2d ) then
00425 !
00426 !===> Store first longitude and latitude corner indices of all donor cells for each target cell 
00427 !
00428          offset = 0
00429          do ipart = 1, search%npart
00430 
00431             do j = 1, npoints(2,ipart)-1
00432                do i = 1, npoints(1,ipart)-1
00433 
00434                   do jj = 1, ic(2,ipart)%vector(j)
00435                      do ii = 1, ic(1,ipart)%vector(i)
00436 
00437                         icell = srclocs(1,ipart)%vector(i) - grid_valid_shape(1,1) + ii 
00438 
00439                         if ( icell <= length(1) ) then
00440                            !
00441                            ! we can do it this way as the corners%vector is already sorted.
00442                            !
00443                            overlap = ( corners(1,ipart)%vector(i+1) - chmin(1)%vector(icell) ) * &
00444                                      ( chmax(1)%vector(icell) - corners(1,ipart)%vector(i) ) >= 0
00445 
00446                            if ( .not. overlap .and. cyclic(1) ) then 
00447 
00448                               icell = srclocs(1,ipart)%vector(i) - grid_valid_shape(1,1) + 1
00449                               if (  corners(1,ipart)%vector(i) < chmin(1)%vector(icell) ) then
00450 #ifdef DEBUGX
00451                                  print *,  trim(ch_id),':',' search cell minimum      ', &
00452                                                              corners(1,ipart)%vector(i) 
00453                                  print *,  trim(ch_id),':',' is outside of local cell ', &
00454                                                              chmin(1)%vector(icell)
00455                                  print *,  trim(ch_id),':',' We look to the left.'
00456 #endif
00457                                  ctr_corner = corners(1,ipart)%vector(i) + global(1)
00458                                  do icell = length(1), 1, -1
00459                                     if ( ctr_corner > chmin(1)%vector(icell) ) exit
00460                                  enddo
00461                                else
00462 #ifdef DEBUGX
00463                                  print *, trim(ch_id),':', ' nothing on the left ', i, j, ii, jj
00464 #endif
00465                                  icell = PSMILe_undef
00466                                endif
00467                            endif
00468                         else
00469 #ifdef DEBUGX
00470                            print *,  trim(ch_id),':',' nothing there ', j, i, jj, ii
00471 #endif
00472                            icell = PSMILe_undef
00473                         endif
00474 
00475                         neighcells_3d(1,offset + ii + (jj-1)*ic(1,ipart)%vector(i), 1 ) &
00476                              = icell
00477                         neighcells_3d(2,offset + ii + (jj-1)*ic(1,ipart)%vector(i), 1 ) &
00478                              = srclocs(2,ipart)%vector(j) + jj - 1
00479 
00480                      enddo
00481                   enddo
00482 
00483                   offset = offset + ic(1,ipart)%vector(i)*ic(2,ipart)%vector(j)
00484 
00485                enddo
00486             enddo
00487 !
00488 !===> Store vertical coordinate indices of all donor points for each target points
00489 !  here for <none> interpolation in the vertical.
00490 !
00491             do k = 1, npoints(3,ipart) ! lenk
00492                neighcells_3d(3,(k-1)*lenij+1:k*lenij,1) = srclocs(3,ipart)%vector(k)
00493             enddo
00494 
00495          enddo !ipart
00496 !
00497 !===> Store horizontal information for all levels k > 1
00498 !
00499          do k = 2, npoints(3,1) ! lenk
00500             neighcells_3d(1,(k-1)*lenij+1:k*lenij,1) = neighcells_3d(1,1:lenij,1)
00501             neighcells_3d(2,(k-1)*lenij+1:k*lenij,1) = neighcells_3d(2,1:lenij,1)
00502          enddo
00503 
00504       else if ( ndim == ndim_3d ) then
00505 !
00506 !===> Store first longitude and latitude corner indices of all donor cells for each target cell 
00507 !
00508          offset = 0
00509 
00510          do ipart = 1, search%npart
00511             do k = 1, npoints(3,ipart)-1
00512                do j = 1, npoints(2,ipart)-1
00513                   do i = 1, npoints(1,ipart)-1
00514 
00515                      do kk = 1, ic(3,ipart)%vector(k)
00516                         do jj = 1, ic(2,ipart)%vector(j)
00517                            do ii = 1, ic(1,ipart)%vector(i)
00518 
00519                               icell = srclocs(1,ipart)%vector(i) + ii - 1
00520 
00521                               if ( icell <= length(1) ) then
00522                                  if ( corners(1,ipart)%vector(i  ) >= chmin(1)%vector(icell) .and. &
00523                                       corners(1,ipart)%vector(i+1) <= chmax(1)%vector(icell) ) then 
00524                                     n = icell
00525                                  endif
00526                               else
00527                                  n = 99
00528                               endif
00529 
00530                               neighcells_3d(1,offset + ii + (jj-1)*ic(1,ipart)%vector(i), 1 ) &
00531                                    = n
00532                               neighcells_3d(2,offset + ii + (jj-1)*ic(1,ipart)%vector(i), 1 ) &
00533                                    = srclocs(2,ipart)%vector(j) + jj - 1
00534                            enddo
00535                         enddo
00536                      enddo
00537 
00538                      offset = offset + ic(1,ipart)%vector(i) &
00539                                      * ic(2,ipart)%vector(j) &
00540                                      * ic(3,ipart)%vector(k)
00541                   enddo
00542                enddo
00543             enddo
00544          enddo
00545 !
00546 !===> Store vertical coordinate indices of all donor points for each target points
00547 !   Expand to store full volume information
00548 !
00549          do k = 1, lenk
00550             neighcells_3d(3,(k-1)*lenij+1:k*lenij,1) = neighcells_3d(3,1:lenij,1)
00551          enddo
00552 
00553       endif ! ndim
00554 !
00555 !===> Store other corner indices for corners 2, 3, and 4.
00556 !
00557       neighcells_3d(1,1:nbr_cells_tot,2) = neighcells_3d(1,1:nbr_cells_tot,1) + 1
00558       neighcells_3d(1,1:nbr_cells_tot,3) = neighcells_3d(1,1:nbr_cells_tot,1) + 1
00559       neighcells_3d(1,1:nbr_cells_tot,4) = neighcells_3d(1,1:nbr_cells_tot,1)
00560 
00561       neighcells_3d(2,1:nbr_cells_tot,2) = neighcells_3d(2,1:nbr_cells_tot,1)
00562       neighcells_3d(2,1:nbr_cells_tot,3) = neighcells_3d(2,1:nbr_cells_tot,1) + 1
00563       neighcells_3d(2,1:nbr_cells_tot,4) = neighcells_3d(2,1:nbr_cells_tot,1) + 1
00564 
00565       neighcells_3d(3,1:nbr_cells_tot,2) = neighcells_3d(3,1:nbr_cells_tot,1)
00566       neighcells_3d(3,1:nbr_cells_tot,3) = neighcells_3d(3,1:nbr_cells_tot,1)
00567       neighcells_3d(3,1:nbr_cells_tot,4) = neighcells_3d(3,1:nbr_cells_tot,1)
00568 !
00569 !===> Store other corner indices for corners 5, 6, 7, and 8 when required, e.g.
00570 !     for 2D1D interpolation with 1D in the vertical.
00571 !
00572       if ( n_corners == 8 ) then
00573          do n = 5, n_corners
00574            neighcells_3d(1,1:nbr_cells_tot,n) = neighcells_3d(1,1:nbr_cells_tot,n-4)
00575            neighcells_3d(2,1:nbr_cells_tot,n) = neighcells_3d(2,1:nbr_cells_tot,n-4)
00576            neighcells_3d(3,1:nbr_cells_tot,n) = neighcells_3d(3,1:nbr_cells_tot,1)+1
00577          enddo
00578       endif
00579 !
00580 !===> Set indices for cyclic coordinate directions
00581 !
00582       do i = 1, ndim_3d
00583 
00584          if (cyclic(i) .and. length(i) > 1) then
00585 
00586             do n = 1, n_corners
00587 
00588                do j = 1, nbr_cells_tot
00589 
00590                   if ( neighcells_3d(i,j,n) < grid_valid_shape(1,i)) then
00591                      neighcells_3d(i,j,n) = neighcells_3d(i,j,n) + length (i)
00592                   else if ( neighcells_3d(i,j,n) > grid_valid_shape(2,i)) then
00593                      neighcells_3d(i,j,n) = neighcells_3d(i,j,n) - length (i)
00594                   endif
00595 
00596                enddo
00597 
00598             end do ! n
00599 
00600          endif
00601 
00602       end do ! j
00603 !
00604 !===> Control special cases (2d hyperplanes)
00605 !
00606       do i = 1, ndim_3d
00607          if (length(i) == 1) then
00608             neighcells_3d(i,1:nbr_cells_tot,:) = grid_valid_shape(1,i)
00609          endif
00610       end do
00611 !
00612 !===> All done
00613 !
00614       do i = 1, ndim
00615          Deallocate( chmin(i)%vector, chmax(i)%vector, STAT=ierror )
00616          if (ierror /= 0) then
00617             ierrp (1) = 2 * length(i)
00618             ierror = PRISM_Error_Dealloc
00619             call psmile_error ( ierror, 'chmin, chmax', &
00620                  ierrp, 1, __FILE__, __LINE__ )
00621             return
00622          endif
00623       enddo
00624 
00625       do ipart = 1, search%npart
00626          do i = 1, ndim
00627             Deallocate(ic(i,ipart)%vector, STAT=ierror )
00628             if (ierror /= 0) then
00629                ierrp (1) = npoints(i,ipart)-1
00630                ierror = PRISM_Error_Dealloc
00631                call psmile_error ( ierror, 'ic vector', &
00632                     ierrp, 1, __FILE__, __LINE__ )
00633                return
00634             endif
00635          enddo
00636       enddo
00637 
00638 !=================================================================================
00639 
00640 #ifdef VERBOSE
00641       print 9980, trim(ch_id)
00642 
00643       call psmile_flushstd
00644 #endif /* VERBOSE */
00645 !
00646 !  Formats:
00647 !
00648 9990 format (1x, a, ': psmile_neigh_cells_3d_reg_real')
00649 9980 format (1x, a, ': psmile_neigh_cells_3d_reg_real: eof')
00650 
00651 #ifdef PRISM_ASSERTION
00652 9970 format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
00653 #endif
00654 
00655       end subroutine psmile_neigh_cells_3d_reg_real