psmile_neigh_cells_3d_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! Copyright 2011, DKRZ, Hamburg, Germany.
00004 ! All rights reserved. Use is subject to OASIS4 license terms.
00005 !-----------------------------------------------------------------------
00006 !BOP
00007 !
00008 #if 0
00009 #define GRIDPOINT 1642
00010 #define META_RANK 2
00011 #define META_POST
00012 #endif
00013 ! #define DEBUGX
00014 #define EXTRA_NORTH_SOUTH_BOUND_CHECK
00015 !
00016 ! !ROUTINE: PSMILe_Neigh_cells_3d_dble
00017 !
00018 ! !INTERFACE:
00019       subroutine psmile_neigh_cells_3d_dble ( use_how,      &
00020                grid_shape, interpolation_mode, cyclic,      &
00021                corner_shape_3d, nbr_corners,                &
00022                corner_x, corner_y, corner_z,                &
00023                search,  control, tgt_cell, tgt_corners,     &
00024                npoints, srclocs, msklocs, ncpl,             &
00025                num_neigh, nbr_cells, ierror )
00026 !
00027 ! !USES:
00028 !
00029       use PRISM_constants
00030       use PSMILe, dummy_interface => PSMILe_Neigh_cells_3d_dble
00031       use psmile_reallocate
00032 #ifdef DEBUG_TRACE
00033       use psmile_debug_trace
00034 #endif
00035 
00036       implicit none
00037 !
00038 ! !INPUT PARAMETERS:
00039 !
00040       Integer, Intent (In)            :: use_how(3)
00041 !
00042 !     container of information  under consideration
00043 !
00044       Integer, Intent (In)            :: grid_shape (2, ndim_3d)
00045 !
00046 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00047 !
00048       Integer, Intent (In)            :: interpolation_mode
00049 !
00050 !     Specifies the interpolation mode with
00051 !        interpolation_mode = 3 : 3D conservative remapping
00052 !                                 (not supported)
00053 !        interpolation_mode = 2 : 2D conservative remapping
00054 !        interpolation_mode = 1 : not supported
00055 !                                 (determines 2 corners)
00056 !
00057       Logical, Intent (In)            :: cyclic (ndim_3d)
00058 !    
00059 !     Does the block have cyclic coordinates in the corresponding direction ?
00060 !
00061       Integer, Intent (In)            :: corner_shape_3d (2,ndim_3d,ndim_3d)
00062 !
00063 !     Shape of local corners -> corner_x, corner_y
00064 !
00065       Integer, Intent (In)            :: nbr_corners(ndim_3d)
00066 !
00067 !     Number of corners of local corners corner_x, corner_y, corner_z
00068 !
00069       Type (Enddef_search), Intent (InOut) :: search
00070 !
00071 !     Info's on coordinates to be searched
00072 !
00073       Integer, Intent (In)            :: control (2, ndim_3d, search%search_data%npart)
00074 !
00075 !     Control range
00076 !
00077       Type (dble_vector), Intent (InOut) :: tgt_cell(ndim_3d)
00078 !
00079 !     Cells for which we search
00080 !
00081       Integer, Intent (In)            :: tgt_corners(ndim_3d)
00082 !
00083 !     Number of corner in the target cell
00084 !
00085       Integer, Intent (In)            :: npoints
00086 !
00087 !     Total number of locations to be transferred
00088 !
00089       Integer, Intent (In)            :: ncpl
00090 !
00091 !     Total number of locations
00092 !
00093       Integer, Intent (In)            :: srclocs (ndim_3d, npoints)
00094 
00095 !     Indices of the (method) grid cell in which the point was found.
00096 
00097       Logical, Intent (In)            :: msklocs (npoints)
00098 
00099 !     Mask of source locations array, only true srclocs need to be considered for target cells.
00100 
00101       Integer, Intent (In)            :: num_neigh
00102 !
00103 !     Last dimension of neighbors corner "neighbors_3d"
00104 !
00105       Double Precision, Intent (In)   :: corner_x (corner_shape_3d (1,1,1):corner_shape_3d(2,1,1), 
00106                                                    corner_shape_3d (1,2,1):corner_shape_3d(2,2,1), 
00107                                                    corner_shape_3d (1,3,1):corner_shape_3d(2,3,1), 
00108                                                    nbr_corners(1))
00109 
00110       Double Precision, Intent (In)   :: corner_y (corner_shape_3d (1,1,2):corner_shape_3d(2,1,2), 
00111                                                    corner_shape_3d (1,2,2):corner_shape_3d(2,2,2), 
00112                                                    corner_shape_3d (1,3,2):corner_shape_3d(2,3,2), 
00113                                                    nbr_corners(2))
00114 
00115       Double Precision, Intent (In)   :: corner_z (corner_shape_3d (1,1,3):corner_shape_3d(2,1,3), 
00116                                                    corner_shape_3d (1,2,3):corner_shape_3d(2,2,3), 
00117                                                    corner_shape_3d (1,3,3):corner_shape_3d(2,3,3), 
00118                                                    nbr_corners(3))
00119 !
00120 !     Local corners
00121 !
00122 !
00123 ! !OUTPUT PARAMETERS:
00124 !
00125       Integer, Intent (Out)           :: nbr_cells(ncpl)
00126 !
00127 !     Number of cells found per target cell
00128 !
00129       Integer, Intent (Out)           :: ierror
00130 
00131 !     Returns the error code of PSMILE_Neigh_cells_irreg2_dble;
00132 !             ierror = 0 : No error
00133 !             ierror > 0 : Severe error
00134 !
00135 ! !DEFINED PARAMETERS:
00136 !
00137 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00138 !               = 2 ** ndim_3d
00139 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00140 !               = 2 ** ndim_2d
00141 !
00142       Integer, Parameter              :: n_corners_3d = 2 * 2 * 2
00143       Integer, Parameter              :: n_corners_2d = 2 * 2
00144       Integer, Parameter              :: n_corners_1d = 2
00145 !
00146 ! !LOCAL VARIABLES
00147 !
00148 !     ... for local cell search
00149 !
00150       Type (line_dble), Allocatable  :: tgt(:), rotated_tgt(:)
00151       Type (line_dble), Allocatable  :: src(:), rotated_src(:)
00152  
00153       Integer, Pointer     :: gauss_stack(:)
00154       Integer, Pointer     :: stack(:,:)
00155       Integer              :: curr_src_cell(2), curr_tgt_cell
00156       Logical              :: is_inital_cell
00157       Logical              :: is_boundary_cell
00158       Integer              :: num_boundary_cells
00159       Integer, Pointer     :: tmp_neighcells_3d(:,:)
00160       Integer, Allocatable :: visited(:,:)
00161       Logical              :: overlap
00162       Logical              :: need_rotation, have_rotated_tgt_cell
00163 
00164       Logical              :: northern_boundary_is_checked
00165       Logical              :: southern_boundary_is_checked
00166 
00167       Integer, Parameter   :: stacksize_inc = 512
00168       Integer, Parameter   :: bnd_size_inc  = 512
00169       Integer              :: nbr_cells_inc
00170 
00171       Integer              :: stacksize
00172 
00173       Integer              :: i, j, n
00174       Integer              :: n_corners
00175       Integer              :: num_found_overlaps
00176 !
00177 !     ... For locations controlled
00178 !
00179       Integer                         :: length (ndim_3d)
00180 !
00181       Integer, Parameter              :: nca (ndim_3d) = (/n_corners_1d, n_corners_2d, n_corners_3d/)
00182 
00183       Integer                         :: grid_id
00184 
00185       Type(Grid), Pointer             :: gp
00186 !
00187 !     ... for error handling
00188 !
00189       Integer, Parameter              :: nerrp = 2
00190       Integer                         :: ierrp (nerrp)
00191 !
00192 #ifdef META_POST
00193       Double Precision                :: meta_shift_x, meta_shift_y
00194       Logical                         :: do_meta_post
00195       Character (len=64)              :: meta_post_string
00196 #endif
00197 #ifdef DEBUGX
00198       double precision, parameter     :: srch_tgt_ubound = 1.0d0
00199       double precision, parameter     :: srch_tgt_lbound = -1.0d0
00200 #endif
00201 !
00202 ! !DESCRIPTION:
00203 !
00204 ! Subroutine "PSMILe_Neigh_cells_3d_dble" searches the cells on the donor
00205 ! grid that overlap with a given target cell.
00206 !
00207 ! Subroutine "PSMILe_Neigh_cells_3d_dble" is designed for (target) grids
00208 ! of type "PRISM_irrlonlatvrt".
00209 !
00210 ! !TO DO:
00211 !
00212 !===> Check what we have to do for linear interpolation (interpolation_mode) in the vertical
00213 !      Is this already covered by a subsequent call to neigh_trili_3d?
00214 !
00215 !===> Cyclic boundary conditions
00216 !
00217 !
00218 ! !REVISION HISTORY:
00219 !
00220 !   Date      Programmer   Description
00221 ! ----------  ----------   -----------
00222 ! 06.06.19    R. Redler    created
00223 ! 24.03.2011  M. Hanke     simplified handling of reduced gauss grids
00224 ! 28.03.2011  M. Hanke     further simplifed this routine
00225 ! 06.04.2011  M. Hanke     added rotation for cells close to the pole
00226 !
00227 !EOP
00228 !----------------------------------------------------------------------
00229 !
00230 !  $Id: psmile_neigh_cells_3d_dble.F90 3182 2011-05-10 07:23:16Z coquart $
00231 !  $Author: coquart $
00232 !
00233    Character(len=len_cvs_string), save :: mycvs = 
00234        '$Id: psmile_neigh_cells_3d_dble.F90 3182 2011-05-10 07:23:16Z coquart $'
00235 !
00236 !----------------------------------------------------------------------
00237 !
00238 !===> Prologue
00239 !
00240 #ifdef VERBOSE
00241       print 9990, trim(ch_id)
00242 
00243       call psmile_flushstd
00244 #endif /* VERBOSE */
00245 !
00246 !===> Control interpolation mode
00247 !
00248       if (interpolation_mode < ndim_2d .or. interpolation_mode > ndim_2d) then
00249          ierrp (1) = interpolation_mode
00250          ierror = PRISM_Error_Internal
00251          call psmile_error ( ierror, 'unsupported interpolation mode', &
00252                              ierrp, 1, __FILE__, __LINE__ )
00253          return
00254       endif
00255 !
00256 !===> Initialization
00257 !
00258       ierror = 0
00259       overlap = .false.
00260       grid_id = use_how(3)
00261       gp => Grids(grid_id)
00262 
00263       shlon(0) =    0.0
00264       shlon(1) = -360.0
00265       shlon(2) =  360.0
00266 
00267       ! set an increment close to ncpl but an
00268       ! integer multiple of some power of two
00269       nbr_cells_inc = ((ncpl/16)+1)*16
00270 
00271       n_corners = nca (interpolation_mode)
00272 
00273       if ( n_corners /= N_CORNERS_2D ) then
00274          ierrp (1) = n_corners
00275          ierror = PRISM_Error_Internal
00276          call psmile_error ( ierror, 'unsupported number of corners', &
00277                              ierrp, 1, __FILE__, __LINE__ )
00278       endif
00279 
00280       length(1) = grid_shape(2,1) - grid_shape(1,1) + 1
00281       length(2) = grid_shape(2,2) - grid_shape(1,2) + 1
00282       length(3) = grid_shape(2,3) - grid_shape(1,3) + 1
00283 
00284 #ifdef PRISM_ASSERTION
00285       ! check whether srclocs only contains valid indices
00286       do n = 1, ncpl
00287          if ( srclocs(1,n) < grid_shape(1,1) .or. &
00288               srclocs(1,n) > grid_shape(2,1) .or. &
00289               srclocs(2,n) < grid_shape(1,2) .or. &
00290               srclocs(2,n) > grid_shape(2,2) .or. &
00291               srclocs(3,n) < grid_shape(1,3) .or. &
00292               srclocs(3,n) > grid_shape(2,3)) then
00293 
00294             print *, 'Incorrect index in n', n, srclocs (:, n)
00295             print *, 'grid shape ', grid_shape
00296             call psmile_assert (__FILE__, __LINE__, &
00297                  'wrong index')
00298          endif
00299       enddo
00300 #endif
00301 
00302 #ifdef META_POST
00303       ! 595 points width for an A4 page, scaled by 180
00304       ! a shift my have to be applied in order to get
00305       ! the boxes on display
00306       if ( global_rank == META_RANK ) then
00307          open ( unit=99, file='cell_search.txt', form='formatted')
00308          write(99,'(a)') 'beginfig(-1)'
00309          write(99,'(a)') 'numeric u;'
00310          write(99,'(a3,f8.4,a1)') 'u:=', 595.0/180, ';'
00311          write(99,'(a)') 'pickup pencircle scaled .1pt;'
00312 
00313          meta_shift_x = 0
00314          meta_shift_y = 0
00315       endif
00316 #endif
00317 
00318       allocate(visited(grid_shape(1,1):grid_shape(2,1),  &
00319                        grid_shape(1,2):grid_shape(2,2)), &
00320                search%boundary_cell(bnd_size_inc,5),     &
00321                tmp_neighcells_3d(3,nbr_cells_inc),       &
00322                tgt(tgt_corners(1)), src(n_corners),      &
00323                rotated_tgt(tgt_corners(1)),              &
00324                rotated_src(n_corners), stat=ierror)
00325 
00326       if ( ierror > 0 ) then
00327          ierrp (1) = ierror
00328          ierrp (2) = length(1) * length(2) + 5*bnd_size_inc + 3*nbr_cells_inc + 8
00329          ierror = PRISM_Error_Alloc
00330          call psmile_error ( ierror, 'various variables', &
00331               ierrp, 2, __FILE__, __LINE__ )
00332          return
00333       endif
00334 
00335       ! initialisation
00336       
00337       visited(:,:)              = 0
00338       search%boundary_cell(:,:) = PSMILe_Undef
00339       nbr_cells(:)              = 0
00340       num_boundary_cells        = 0
00341       num_found_overlaps        = 0
00342       curr_tgt_cell             = 0
00343 
00344       ! if the local grid is a reduced gauss grid
00345       if ( gp%grid_type == PRISM_Gaussreduced_regvrt ) then
00346 
00347          allocate(gauss_stack(stacksize_inc), stat=ierror)
00348 
00349          if ( ierror > 0 ) then
00350             ierrp (1) = ierror
00351             ierrp (2) = stacksize_inc
00352             ierror = PRISM_Error_Alloc
00353             call psmile_error ( ierror, 'gauss_stack', &
00354                  ierrp, 2, __FILE__, __LINE__ )
00355             return
00356          endif
00357 
00358          ! for all target points
00359          do n = 1, npoints
00360 
00361             ! msklocs contains source cells that shall not be considered
00362             ! as starting cells for further local search. In particular
00363             ! this is required for reglonlatvrt source cell grids.
00364 
00365             if ( .not. msklocs(n) ) cycle
00366 
00367             curr_tgt_cell  = curr_tgt_cell + 1
00368 
00369 #ifdef META_POST
00370             do_meta_post = curr_tgt_cell == GRIDPOINT .and. global_rank == META_RANK
00371 #endif
00372 
00373             call get_tgt_cell (tgt, tgt_corners(1), curr_tgt_cell)
00374 
00375             have_rotated_tgt_cell = .false.
00376             is_boundary_cell = .false.
00377 
00378             ! add the initial cell to the stack
00379             is_inital_cell         = .true.
00380             stacksize              = 0
00381 
00382             call add_gauss_neigh_to_stack (srclocs(1,n), stacksize)
00383 
00384             stacksize              = stacksize + 1
00385             gauss_stack(stacksize) = srclocs(1,n)
00386 
00387             ! as long as there are cells on the stack
00388             do while (stacksize > 0)
00389 
00390                ! get the top element from the stack
00391                curr_src_cell(1) = gauss_stack(stacksize)
00392                stacksize        = stacksize - 1
00393 
00394                ! if the current source cell is a remote cell
00395                if (curr_src_cell(1) < grid_shape(1,1) .or. &
00396                    curr_src_cell(1) > grid_shape(2,1) ) then
00397 
00398                   is_boundary_cell = .true.
00399                   cycle
00400                endif
00401 
00402               ! if the cell was already visited
00403                if ( visited(curr_src_cell(1),1) == curr_tgt_cell ) then
00404                   cycle
00405                else
00406                   ! mark cell as being visited
00407                   visited(curr_src_cell(1),1) = curr_tgt_cell
00408                endif
00409 
00410                ! get the coordinates of the current source cell
00411                call get_src_cell (src, rotated_src, n_corners, curr_src_cell, &
00412                                   need_rotation)
00413 
00414                ! if the source cell was rotated because it was close to the pole
00415                if (need_rotation) then
00416 
00417                   ! do we already have a rotated version of the current tgt cell?
00418                   if (.not. have_rotated_tgt_cell) then
00419 
00420                      call rotate_tgt(tgt, rotated_tgt, tgt_corners(1), rotated_src, &
00421                                      n_corners)
00422                      have_rotated_tgt_cell = .true.
00423                   endif
00424 
00425                   call psmile_overlap_dble ( n_corners, tgt_corners(1), rotated_src, &
00426                                              rotated_tgt, overlap )
00427                else
00428                   call psmile_overlap_dble ( n_corners, tgt_corners(1), src, tgt, overlap )
00429                endif
00430 
00431 #ifdef META_POST
00432                if ( do_meta_post ) then
00433                   
00434                   if (is_inital_cell) then
00435                      meta_shift_x = -tgt(1)%p1%x+180.0/6.0
00436                      meta_shift_y = -tgt(1)%p1%y+180.0/6.0
00437                   endif
00438 
00439                   if (overlap) then
00440                      meta_post_string = ' source cell point o'
00441                   else
00442                      meta_post_string = ' source cell point v'
00443                   endif
00444                   call mp_fill  ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1, overlap )
00445 
00446                   do j = 1, n_corners
00447                      print *, trim(meta_post_string), curr_src_cell(1), curr_tgt_cell, &
00448                               src(j)%p1%x, src(j)%p1%y
00449                      call mp_draw ( src(j)%p1, src(j)%p2, .false. )
00450                   enddo
00451                   print *, ' --------- '
00452                endif
00453 #endif
00454                ! if there was an overlap
00455                if (overlap) then
00456 
00457                   if (num_found_overlaps+1 > size(tmp_neighcells_3d,2)) then
00458                      tmp_neighcells_3d => psmile_realloc(tmp_neighcells_3d, &
00459                                                          (/3,num_found_overlaps+nbr_cells_inc/))
00460                   endif
00461 
00462                   num_found_overlaps = num_found_overlaps + 1
00463 
00464                   ! save the cell in neighcells_3d
00465                   tmp_neighcells_3d(1,num_found_overlaps) = curr_src_cell(1)
00466                   tmp_neighcells_3d(2,num_found_overlaps) = 1
00467                   tmp_neighcells_3d(3,num_found_overlaps) = srclocs(3,n)
00468 
00469                   nbr_cells(curr_tgt_cell) = nbr_cells(curr_tgt_cell) + 1
00470                endif
00471 
00472                ! if there was an overlap or if this not is the first cell
00473                if (overlap .and. .not. is_inital_cell) then
00474 
00475                   ! add the neighbours of the current cell to the stack
00476                   call add_gauss_neigh_to_stack (curr_src_cell(1), stacksize)
00477                endif
00478 
00479                is_inital_cell = .false.
00480             enddo ! (stacksize > 0)
00481 
00482             ! if the current target cell is at the boundary of the local source grid
00483             if ( is_boundary_cell .and. nbr_cells(curr_tgt_cell) > 0) then
00484 
00485                if (num_boundary_cells + 1 > size(search%boundary_cell,1)) then
00486 
00487                   search%boundary_cell => psmile_realloc(search%boundary_cell, &
00488                                                          (/size(search%boundary_cell,1) + &
00489                                                          bnd_size_inc,5/))
00490                   search%boundary_cell(num_boundary_cells+1:,:) = PSMILe_Undef
00491                endif
00492 
00493                num_boundary_cells = num_boundary_cells + 1
00494 
00495                search%boundary_cell(num_boundary_cells,1) = curr_tgt_cell
00496                search%boundary_cell(num_boundary_cells,5) = psmile_rank
00497             endif
00498 
00499 #ifdef META_POST
00500             if ( do_meta_post ) then
00501                do j = 1, tgt_corners(1)
00502                   call mp_draw ( tgt(j)%p1, tgt(j)%p2, .true. )
00503                   print *, ' target cell point ', curr_tgt_cell, &
00504                            tgt(j)%p1%x, tgt(j)%p1%y
00505                enddo
00506                print *, ' --------- '
00507             endif
00508 #endif
00509          enddo ! n loop over npoints
00510 
00511          deallocate (gauss_stack)
00512 
00513       else ! PRISM_Gaussreduced_regvrt
00514 
00515          allocate(stack(2,stacksize_inc), stat=ierror)
00516 
00517          if ( ierror > 0 ) then
00518             ierrp (1) = ierror
00519             ierrp (2) = nbr_cells_inc
00520             ierror = PRISM_Error_Alloc
00521             call psmile_error ( ierror, 'stack', ierrp, 2, &
00522                                 __FILE__, __LINE__ )
00523             return
00524          endif
00525 
00526          ! for all target points
00527          do n = 1, npoints
00528 
00529             ! msklocs contains source cells that shall not be considered
00530             ! as starting cells for further local search. In particular
00531             ! this is required for reglonlatvrt source cell grids.
00532 
00533             if ( .not. msklocs(n) ) cycle
00534 
00535             curr_tgt_cell  = curr_tgt_cell + 1
00536 
00537 #ifdef META_POST
00538             do_meta_post = curr_tgt_cell == GRIDPOINT .and. global_rank == META_RANK
00539 #endif
00540             call get_tgt_cell (tgt, tgt_corners(1), curr_tgt_cell)
00541 
00542             have_rotated_tgt_cell = .false.
00543             is_boundary_cell = .false.
00544 
00545             ! add the initial cell to the stack
00546             is_inital_cell = .true.
00547             stacksize      = 0
00548 
00549             call add_neigh_to_stack (srclocs(1,n), srclocs(2,n), stacksize)
00550 
00551             stacksize      = stacksize + 1
00552 
00553             ! adjust size of stack if necessary
00554             if (size(stack,2) < stacksize) then
00555                stack => psmile_realloc(stack, (/2,size(stack) + stacksize_inc/))
00556             endif
00557 
00558             stack(1, stacksize) = srclocs(1,n)
00559             stack(2, stacksize) = srclocs(2,n)
00560 
00561 #ifdef EXTRA_NORTH_SOUTH_BOUND_CHECK
00562             northern_boundary_is_checked = .false.
00563             southern_boundary_is_checked = .false.
00564 #else
00565             northern_boundary_is_checked = .true.
00566             southern_boundary_is_checked = .true.
00567 #endif
00568 
00569             ! as long as there are cells on the stack
00570             do while (stacksize > 0)
00571 
00572                ! get the top element from the stack
00573                curr_src_cell(:) = stack(:,stacksize)
00574                stacksize        = stacksize - 1
00575 
00576                ! if the current source cell is a remote cell
00577                if ( curr_src_cell(1) < grid_shape(1,1) .or. &
00578                     curr_src_cell(1) > grid_shape(2,1) .or. &
00579                     curr_src_cell(2) < grid_shape(1,2) .or. &
00580                     curr_src_cell(2) > grid_shape(2,2) ) then
00581 
00582                   if ( gp%grid_type == PRISM_Irrlonlat_regvrt ) then
00583                   
00584                      if ( curr_src_cell(2) < grid_shape(1,2) .and. &
00585                          .not. southern_boundary_is_checked) then
00586 
00587                         southern_boundary_is_checked = .true.
00588 
00589                         call add_boundary_to_stack(.false., stacksize)
00590                      endif
00591 
00592                      if ( curr_src_cell(2) > grid_shape(2,2) .and. &
00593                           .not. northern_boundary_is_checked) then
00594 
00595                         northern_boundary_is_checked = .true.
00596 
00597                         call add_boundary_to_stack(.true., stacksize)
00598                      endif
00599                   endif
00600 
00601                   is_boundary_cell = .true.
00602                   cycle
00603                endif
00604 
00605               ! if the cell was already visited
00606                if ( visited(curr_src_cell(1),curr_src_cell(2)) == curr_tgt_cell ) then
00607                   cycle
00608                else
00609                   ! mark cell as being visited
00610                   visited(curr_src_cell(1),curr_src_cell(2)) = curr_tgt_cell
00611                endif
00612 
00613                ! get the coordinates of the current source cell
00614                call get_src_cell(src, rotated_src, n_corners, curr_src_cell, need_rotation)
00615 
00616                ! if the source cell was rotated because it was close to the pole
00617                if (need_rotation) then
00618 
00619                   ! do we already have a rotated version of the current tgt cell?
00620                   if (.not. have_rotated_tgt_cell) then
00621 
00622                      call rotate_tgt(tgt, rotated_tgt, tgt_corners(1), rotated_src, &
00623                                      n_corners)
00624                      have_rotated_tgt_cell = .true.
00625                   endif
00626 
00627                   call psmile_overlap_dble ( n_corners, tgt_corners(1), rotated_src, &
00628                                              rotated_tgt, overlap )
00629                else
00630                   call psmile_overlap_dble ( n_corners, tgt_corners(1), src, tgt, overlap )
00631                endif
00632 
00633 #ifdef META_POST
00634                if ( do_meta_post ) then
00635                   
00636                   if (is_inital_cell) then
00637                      meta_shift_x = -tgt(1)%p1%x+180.0/6.0
00638                      meta_shift_y = -tgt(1)%p1%y+180.0/6.0
00639                   endif
00640 
00641                   if (overlap) then
00642                      meta_post_string = ' source cell point o'
00643                   else
00644                      meta_post_string = ' source cell point v'
00645                   endif
00646                   call mp_fill  ( src(1)%p1, src(2)%p1, src(3)%p1, src(4)%p1, overlap )
00647 
00648                   do j = 1, n_corners
00649                      print *, trim(meta_post_string), curr_src_cell(:), curr_tgt_cell, &
00650                               src(j)%p1%x, src(j)%p1%y
00651                      call mp_draw ( src(j)%p1, src(j)%p2, .false. )
00652                   enddo
00653                   print *, ' --------- '
00654                endif
00655 #endif
00656                ! if there was an overlap
00657                if (overlap) then
00658 
00659                   if (num_found_overlaps+1 > size(tmp_neighcells_3d,2)) then
00660                      tmp_neighcells_3d => psmile_realloc(tmp_neighcells_3d, &
00661                                                          (/3,num_found_overlaps+nbr_cells_inc/))
00662                   endif
00663 
00664                   num_found_overlaps = num_found_overlaps + 1
00665 
00666                   ! save the cell in neighcells_3d
00667                   tmp_neighcells_3d(1:2,num_found_overlaps) = curr_src_cell(1:2)
00668                   tmp_neighcells_3d(3,num_found_overlaps)   = srclocs(3,n)
00669 
00670                   nbr_cells(curr_tgt_cell)                  = nbr_cells(curr_tgt_cell) + 1
00671                endif
00672 
00673                ! if there was an overlap or if this not is the first cell
00674                if (overlap .and. .not. is_inital_cell) then
00675 
00676                   ! add the neighbours of the current cell to the stack
00677                   call add_neigh_to_stack (curr_src_cell(1), curr_src_cell(2), stacksize)
00678                endif
00679 
00680                is_inital_cell = .false.
00681             enddo ! (stacksize > 0)
00682 
00683             ! if the current target cell is at the boundary of the local source grid
00684             if ( is_boundary_cell .and. nbr_cells(curr_tgt_cell) > 0) then
00685 
00686                if (num_boundary_cells + 1 > size(search%boundary_cell,1)) then
00687 
00688                   search%boundary_cell => psmile_realloc(search%boundary_cell, &
00689                                                          (/size(search%boundary_cell,1) + &
00690                                                          bnd_size_inc,5/))
00691 
00692                   search%boundary_cell(num_boundary_cells+1:,:) = PSMILe_Undef
00693                endif
00694 
00695                num_boundary_cells = num_boundary_cells + 1
00696 
00697                search%boundary_cell(num_boundary_cells,1) = curr_tgt_cell
00698                search%boundary_cell(num_boundary_cells,5) = psmile_rank
00699             endif
00700 
00701 #ifdef META_POST
00702             if ( do_meta_post ) then
00703                do j = 1, tgt_corners(1)
00704                   call mp_draw ( tgt(j)%p1, tgt(j)%p2, .true. )
00705                   print *, ' target cell point ', curr_tgt_cell, &
00706                            tgt(j)%p1%x, tgt(j)%p1%y
00707                enddo
00708                print *, ' --------- '
00709             endif
00710 #endif
00711          enddo ! n loop over npoints
00712 
00713          deallocate (stack)
00714 
00715       endif ! Gaussreduced_regvrt
00716 
00717 #ifdef META_POST
00718       if ( global_rank == META_RANK ) then
00719          write(99,'(a)') 'endfig'
00720          write(99,'(a)') 'end'
00721          close(unit=99)
00722       endif
00723 #endif
00724       !
00725       ! Copy temporary results to final destination and deallocate temporary buffer
00726       !
00727       select case ( gp%grid_type )
00728 
00729       case ( PRISM_Reglonlatvrt )
00730 
00731          Allocate ( neighcells_3d(ndim_3d,num_found_overlaps,n_corners), STAT=ierror)
00732 
00733          if (ierror /= 0) then
00734             ierrp (1) = 3 * num_found_overlaps * n_corners
00735             ierror = PRISM_Error_Alloc
00736             call psmile_error ( ierror, 'neighcells_3d', &
00737                  ierrp, 1, __FILE__, __LINE__ )
00738             return
00739          endif
00740 
00741          neighcells_3d(1:3,1:num_found_overlaps,1) = tmp_neighcells_3d(1:3,1:num_found_overlaps)
00742 
00743          neighcells_3d(1,1:num_found_overlaps,2) = neighcells_3d(1,1:num_found_overlaps,1) + 1
00744          neighcells_3d(1,1:num_found_overlaps,3) = neighcells_3d(1,1:num_found_overlaps,1) + 1
00745          neighcells_3d(1,1:num_found_overlaps,4) = neighcells_3d(1,1:num_found_overlaps,1)
00746 
00747          neighcells_3d(2,1:num_found_overlaps,2) = neighcells_3d(2,1:num_found_overlaps,1)
00748          neighcells_3d(2,1:num_found_overlaps,3) = neighcells_3d(2,1:num_found_overlaps,1) + 1
00749          neighcells_3d(2,1:num_found_overlaps,4) = neighcells_3d(2,1:num_found_overlaps,1) + 1
00750 
00751          neighcells_3d(3,1:num_found_overlaps,2) = neighcells_3d(3,1:num_found_overlaps,1)
00752          neighcells_3d(3,1:num_found_overlaps,3) = neighcells_3d(3,1:num_found_overlaps,1)
00753          neighcells_3d(3,1:num_found_overlaps,4) = neighcells_3d(3,1:num_found_overlaps,1)
00754          !
00755          !===> Store other corner indices for corners 5, 6, 7, and 8 when required, e.g.
00756          !     for 2D1D interpolation with 1D in the vertical.
00757          !
00758          if ( n_corners == 8 ) then
00759             do n = 5, n_corners
00760                neighcells_3d(1,1:num_found_overlaps,n) = neighcells_3d(1,1:num_found_overlaps,n-4)
00761                neighcells_3d(2,1:num_found_overlaps,n) = neighcells_3d(2,1:num_found_overlaps,n-4)
00762                neighcells_3d(3,1:num_found_overlaps,n) = neighcells_3d(3,1:num_found_overlaps,1)+1
00763             enddo
00764          endif
00765          !
00766          !===> Set indices for cyclic coordinate directions
00767          !
00768          do i = 1, ndim_3d
00769 
00770             if (cyclic(i) .and. length(i) > 1) then
00771 
00772                do n = 1, n_corners
00773 
00774                   do j = 1, num_found_overlaps
00775 
00776                      if ( neighcells_3d(i,j,n) < grid_shape(1,i)) then
00777                           neighcells_3d(i,j,n) = neighcells_3d(i,j,n) + length (i)
00778                      else if ( neighcells_3d(i,j,n) > grid_shape(2,i)) then
00779                           neighcells_3d(i,j,n) = neighcells_3d(i,j,n) - length (i)
00780                      endif
00781 
00782                   enddo
00783 
00784                end do ! n
00785 
00786             endif
00787 
00788          end do ! j
00789          !
00790          !===> Control special cases (2d hyperplanes)
00791          !
00792          do i = 1, ndim_3d
00793             if (length(i) == 1) then
00794                neighcells_3d(i,1:num_found_overlaps,:) = grid_shape(1,i)
00795             endif
00796          end do
00797 
00798       case DEFAULT
00799 
00800          allocate(neighcells_3d(3,num_found_overlaps,1), STAT=ierror)
00801 
00802          if ( ierror > 0 ) then
00803             ierrp (1) = ierror
00804             ierrp (2) = 3 * num_found_overlaps
00805             ierror = PRISM_Error_Alloc
00806             call psmile_error ( ierror, 'neighcells_3d', &
00807                  ierrp, 2, __FILE__, __LINE__ )
00808             return
00809          endif
00810 
00811          neighcells_3d(1:3,1:num_found_overlaps,1) = tmp_neighcells_3d(1:3,1:num_found_overlaps)
00812 
00813       end select ! ( gp%grid_type )
00814 
00815       if ( num_boundary_cells > 0 ) then
00816 
00817          search%boundary_cell => psmile_realloc(search%boundary_cell, (/num_boundary_cells,5/))
00818 
00819          ! We should memorise a cell only when neighbours have been found
00820 
00821          do i = 1, num_boundary_cells
00822             if ( search%boundary_cell(i,1) /= PSMILe_Undef ) then
00823                if (nbr_cells(search%boundary_cell(i,1)) == 0 ) &
00824                     search%boundary_cell(i,1) = PSMILe_Undef
00825             endif
00826 #ifdef DEBUGX
00827             !
00828             ! boundary cell indices 2, 3 and 4 are set in psmile_info_trs_loc_*
00829             !
00830             if ( search%boundary_cell(i,1) /= PSMILe_Undef ) then
00831                 print '(a,4i8)', 'boundary cell info ', &
00832                    i, search%boundary_cell(i,1), search%boundary_cell(i,5), &
00833                    nbr_cells(search%boundary_cell(i,1))
00834             endif
00835 #endif /* DEBUGX */
00836          enddo
00837       else
00838 
00839          deallocate(search%boundary_cell, STAT=ierror)
00840 
00841          if ( ierror > 0 ) then
00842             ierrp (1) = ierror
00843             ierrp (2) = 3 * num_found_overlaps
00844             ierror = PRISM_Error_Dealloc
00845             call psmile_error ( ierror, 'search%boundary_cell', ierrp, 2, &
00846                                 __FILE__, __LINE__ )
00847             return
00848          endif
00849 
00850       endif ! ( num_boundary_cells > 0 )
00851 !
00852 !===> Deallocate local memory
00853 !
00854       deallocate(visited, tgt, src, rotated_tgt, rotated_src, &
00855                  tmp_neighcells_3d, stat=ierror)
00856 
00857       if ( ierror > 0 ) then
00858          ierrp (1) = ierror
00859          ierrp (2) = product(length(1:2)) + tgt_corners(1) + n_corners + 3 * num_found_overlaps
00860          ierror = PRISM_Error_Dealloc
00861          call psmile_error ( ierror, 'visited, tmp_neighcells_3d, tgt and src', ierrp, 2, &
00862                              __FILE__, __LINE__ )
00863          return
00864       endif
00865 
00866 #ifdef DEBUG_TRACE
00867       if (ictl > 0 .and. ictl <= ncpl) then
00868          print *, "ictl nbr_cells(ictl)", nbr_cells(ictl)
00869          do i = 0, nbr_cells(ictl)-1
00870             print *, "ictl local neighcell index", &
00871                      neighcells_3d(:,sum(nbr_cells(1:ictl))+i,:)
00872          enddo
00873       endif
00874 #endif
00875 !
00876 !===> All done
00877 !
00878 #ifdef VERBOSE
00879       print 9980, trim(ch_id)
00880 
00881       call psmile_flushstd
00882 #endif /* VERBOSE */
00883 !
00884 !  Formats:
00885 !
00886 9990 format (1x, a, ': psmile_neigh_cells_3d_dble')
00887 9980 format (1x, a, ': psmile_neigh_cells_3d_dble: eof')
00888 
00889       contains
00890 
00891          subroutine get_tgt_cell (tgt, num_corners, tgt_cell_index)
00892 
00893             implicit none
00894 
00895             integer, intent (in)              :: num_corners, tgt_cell_index
00896             type (line_dble), intent (inout)  :: tgt(num_corners)
00897 
00898             integer          ::  i, i_
00899             double precision :: tgt_max, tgt_min
00900 
00901             do i = 0, num_corners-1
00902                i_ = mod(i+1,num_corners)
00903                tgt(i+1)%p1%x = tgt_cell(1)%vector(i *ncpl+tgt_cell_index)
00904                tgt(i+1)%p1%y = tgt_cell(2)%vector(i *ncpl+tgt_cell_index)
00905                tgt(i+1)%p2%x = tgt_cell(1)%vector(i_*ncpl+tgt_cell_index)
00906                tgt(i+1)%p2%y = tgt_cell(2)%vector(i_*ncpl+tgt_cell_index)
00907             enddo
00908             !
00909             ! convert target longitudes to [0,360] interval if necessary
00910             !
00911             tgt_min =  999.0d0
00912             tgt_max = -999.0d0
00913 
00914             do i = 1, num_corners
00915               tgt_min = min(tgt_min,tgt(i)%p1%x)
00916               tgt_max = max(tgt_max,tgt(i)%p1%x)
00917             enddo
00918 
00919 #ifdef DEBUGX
00920             !
00921             ! Find a particular entry
00922             !
00923             if ( tgt(1)%p1%y > srch_tgt_lbound .and. &
00924                  tgt(1)%p1%y < srch_tgt_ubound ) then
00925                print *, ' selected x Corners for ', tgt_cell_index, tgt(:)%p1%x
00926                print *, ' selected y Corners for ', tgt_cell_index, tgt(:)%p1%y
00927             endif
00928 #endif /* DEBUGX */
00929 
00930             if ( tgt_max - tgt_min > 180.0d0 ) then
00931 
00932                do i = 1, num_corners
00933 
00934                   if ( tgt(i)%p1%x > 360.0d0 ) then
00935                        tgt(i)%p1%x = tgt(i)%p1%x - 360.0d0
00936                        tgt_cell(1)%vector((i-1)*ncpl+tgt_cell_index) = &
00937                        tgt_cell(1)%vector((i-1)*ncpl+tgt_cell_index) - 360.0d0
00938                   endif
00939 
00940                   if ( tgt(i)%p1%x < zero    ) then
00941                        tgt(i)%p1%x = tgt(i)%p1%x + 360.0d0
00942                        tgt_cell(1)%vector((i-1)*ncpl+tgt_cell_index) = &
00943                        tgt_cell(1)%vector((i-1)*ncpl+tgt_cell_index) + 360.0d0
00944                   endif
00945 
00946                   i_ = mod(i,num_corners)
00947 
00948                   if ( tgt(i)%p2%x > 360.0d0 ) then
00949                        tgt(i)%p2%x = tgt(i)%p2%x - 360.0d0
00950                        tgt_cell(1)%vector(i_*ncpl+tgt_cell_index) = &
00951                        tgt_cell(1)%vector(i_*ncpl+tgt_cell_index) - 360.0d0
00952                   endif
00953 
00954                   if ( tgt(i)%p2%x < zero    ) then
00955                       tgt(i)%p2%x = tgt(i)%p2%x + 360.0d0
00956                       tgt_cell(1)%vector(i_*ncpl+tgt_cell_index) = &
00957                       tgt_cell(1)%vector(i_*ncpl+tgt_cell_index) + 360.0d0
00958                   endif
00959 
00960                enddo
00961             endif
00962          end subroutine get_tgt_cell
00963 
00964 !-------------------------------------------------------------------------------------
00965 
00966          subroutine get_src_cell (src, rotated_src, num_corners, src_cell_index, &
00967                                   need_rotation)
00968 
00969             use psmile_grid, only : common_grid_range, pole_threshold, &
00970                                     psmile_transform_cell_cyclic
00971 
00972             implicit none
00973             
00974             integer, intent (in)           :: num_corners, src_cell_index(2)
00975             type (line_dble), intent (out) :: src(num_corners), rotated_src(num_corners)
00976             logical, intent(out)           :: need_rotation
00977 
00978             integer            :: j, j_
00979             double precision   :: period
00980             integer, parameter :: edge(2,4) = reshape ((/1,1,2,1,2,2,1,2/), (/2,4/))
00981             need_rotation = .false.
00982 
00983             if ( gp%grid_type == PRISM_Reglonlatvrt       .or. &
00984                  gp%grid_type == PRISM_Reglonlat_sigmavrt) then
00985 
00986                do j = 1, num_corners
00987                   j_ = mod(j,num_corners)+1
00988                   src(j)%p1%x = corner_x(src_cell_index(1),               1,1,edge(1,j))
00989                   src(j)%p1%y = corner_y(               1,src_cell_index(2),1,edge(2,j))
00990                   src(j)%p2%x = corner_x(src_cell_index(1),               1,1,edge(1,j_))
00991                   src(j)%p2%y = corner_y(               1,src_cell_index(2),1,edge(2,j_))
00992                enddo
00993             else if (gp%grid_type == PRISM_Gaussreduced_regvrt) then
00994 
00995                do j = 1, num_corners
00996                   j_ = mod(j,num_corners)+1
00997                   src(j)%p1%x = corner_x(src_cell_index(1),               1,1,edge(1,j))
00998                   src(j)%p1%y = corner_y(               1,src_cell_index(1),1,edge(2,j))
00999                   src(j)%p2%x = corner_x(src_cell_index(1),               1,1,edge(1,j_))
01000                   src(j)%p2%y = corner_y(               1,src_cell_index(1),1,edge(2,j_))
01001                enddo
01002             else
01003                do j = 1, num_corners
01004                   j_ = mod(j,num_corners)+1
01005                   src(j)%p1%x = corner_x(src_cell_index(1),src_cell_index(2),1,j)
01006                   src(j)%p1%y = corner_y(src_cell_index(1),src_cell_index(2),1,j)
01007                   src(j)%p2%x = corner_x(src_cell_index(1),src_cell_index(2),1,j_)
01008                   src(j)%p2%y = corner_y(src_cell_index(1),src_cell_index(2),1,j_)
01009                enddo
01010             endif
01011 
01012             ! if we are close to the pole we have to rotate the cell in order to be
01013             ! able to handle distorted cell that are very close to or on the pole
01014             ! (this is no issue for reglonlat because horizontal edges are always on
01015             ! latitude circles and verticals ones on great circles throuth the poles)
01016             if (minval(src(:)%p1%y) > pole_threshold .or. &
01017                 maxval(src(:)%p1%y) < -pole_threshold) then
01018 
01019                need_rotation = .true.
01020 
01021                ! translate source grid coordintates to the equator
01022                do j = 1, num_corners
01023 
01024                   call psmile_transrot_dble (        src(j)%p1%x,         src(j)%p1%y, &
01025                                              rotated_src(j)%p1%x, rotated_src(j)%p1%y)
01026                enddo ! j
01027 
01028                ! the transformation might cause some bad coordinates
01029                period = common_grid_range(2,1) - common_grid_range(1,1)
01030                call psmile_transform_cell_cyclic ( rotated_src(:)%p1%x, period, ierror )
01031 
01032                do j = 1, num_corners
01033 
01034                   j_ = mod(j,num_corners)+1
01035                   rotated_src(j)%p2 = rotated_src(j_)%p1
01036                enddo ! j
01037             endif
01038          end subroutine get_src_cell
01039 
01040 !-------------------------------------------------------------------------------------
01041 
01042          subroutine rotate_tgt (tgt, rotated_tgt, tgt_num_corners, &
01043                                 rotated_src, src_num_corners)
01044 
01045             use psmile_grid, only : common_grid_range, psmile_transform_cell_cyclic
01046 
01047             implicit none
01048 
01049             integer, intent (in)           :: tgt_num_corners, src_num_corners
01050             type (line_dble), intent (in)  :: rotated_src(src_num_corners), 
01051                                               tgt(tgt_num_corners)
01052             type (line_dble), intent (out) :: rotated_tgt(tgt_num_corners)
01053 
01054             integer          :: j, j_, ierror
01055             double precision :: period, period2
01056 
01057             period = common_grid_range(2,1) - common_grid_range(1,1)
01058             period2 = period * 0.5d0
01059 
01060             ! translate target grid coordintates to the equator
01061             do j = 1, tgt_num_corners
01062 
01063                call psmile_transrot_dble (        tgt(j)%p1%x,         tgt(j)%p1%y, &
01064                                           rotated_tgt(j)%p1%x, rotated_tgt(j)%p1%y)
01065             enddo ! j
01066 
01067             ! the transformation might cause some bad coordinates
01068             call psmile_transform_cell_cyclic ( rotated_tgt(:)%p1%x, period, ierror )
01069          
01070             ! move target longitudes into the common local longitude range
01071             do while ( minval(rotated_tgt(:)%p1%x) < minval(rotated_src(:)%p1%x) - period2 )
01072                rotated_tgt(:)%p1%x =  rotated_tgt(:)%p1%x + period
01073             enddo
01074          
01075             do while ( maxval(rotated_tgt(:)%p1%x) > maxval(rotated_src(:)%p1%x) + period2 )
01076                rotated_tgt(:)%p1%x =  rotated_tgt(:)%p1%x - period
01077             enddo
01078 
01079             do j = 1, tgt_num_corners
01080 
01081                j_ = mod(j,tgt_num_corners)+1
01082                rotated_tgt(j)%p2 = rotated_tgt(j_)%p1
01083             enddo ! j
01084          end subroutine rotate_tgt
01085 
01086 !-------------------------------------------------------------------------------------
01087 
01088          subroutine add_gauss_neigh_to_stack (base, stack_entries)
01089 
01090             use psmile_grid_reduced_gauss
01091 
01092             implicit none
01093 
01094             integer, intent(in)             :: base
01095             integer, intent(inout)          :: stack_entries
01096 
01097             integer          :: base_3d(ndim_3d)
01098             integer, pointer :: stencil_lookup(:)
01099             integer, pointer :: nbr_points_per_lat(:)
01100 
01101             stencil_lookup => gp%reduced_gauss_data%local_1d_stencil_lookup(:,base)
01102             nbr_points_per_lat => gp%reduced_gauss_data%nbr_points_per_lat
01103             
01104             ! adjust size of stack if necessary
01105             if (size(gauss_stack) < stack_entries + 8) then
01106                gauss_stack => psmile_realloc(gauss_stack, size(gauss_stack) + stacksize_inc)
01107             endif
01108             
01109             ! get the global 3d index of the base
01110             base_3d = psmile_gauss_local_1d_to_3d (grid_id, base)
01111 
01112             ! add the eastern and western neighbours to the stack
01113             gauss_stack(stack_entries+1) = stencil_lookup(5)
01114             gauss_stack(stack_entries+2) = stencil_lookup(7)
01115             stack_entries = stack_entries + 2
01116 
01117             ! add the northern neighbour(s)
01118             if (nbr_points_per_lat(base_3d(2)) /= &
01119                nbr_points_per_lat(min(base_3d(2)+1,size(nbr_points_per_lat)))) then
01120 
01121                gauss_stack(stack_entries+1) = stencil_lookup(9)
01122                gauss_stack(stack_entries+2) = stencil_lookup(10)
01123                gauss_stack(stack_entries+3) = stencil_lookup(11)
01124                stack_entries = stack_entries + 3
01125             else
01126                gauss_stack(stack_entries+1) = stencil_lookup(10)
01127                stack_entries = stack_entries + 1
01128             endif
01129 
01130             ! add the southern neighbour(s)
01131             if (nbr_points_per_lat(base_3d(2)) /= &
01132                nbr_points_per_lat(max(base_3d(2)-1,1))) then
01133 
01134                gauss_stack(stack_entries+1) = stencil_lookup(1)
01135                gauss_stack(stack_entries+2) = stencil_lookup(2)
01136                gauss_stack(stack_entries+3) = stencil_lookup(3)
01137                stack_entries = stack_entries + 3
01138             else
01139                gauss_stack(stack_entries+1) = stencil_lookup(2)
01140                stack_entries = stack_entries + 1
01141             endif
01142          end subroutine add_gauss_neigh_to_stack
01143 
01144 !-------------------------------------------------------------------------------------
01145 
01146          subroutine add_neigh_to_stack (base_i, base_j, stack_entries)
01147 
01148             implicit none
01149 
01150             integer, intent(in)             :: base_i, base_j
01151             integer, intent(inout)          :: stack_entries
01152 
01153             integer :: neighbors_i(4)
01154             integer :: neighbors_j(4)
01155 
01156             ! adjust size of stack if necessary
01157             if (size(stack, 2) < stack_entries + 4) then
01158                stack => psmile_realloc(stack, (/2,size(stack) + stacksize_inc/))
01159             endif
01160 
01161             ! compute the indices of the four neighbours
01162             neighbors_i(:) = base_i + (/-1, 0,1,0/)
01163             neighbors_j(:) = base_j + (/ 0,-1,0,1/)
01164 
01165             ! apply cyclic boundary conditions if required
01166             if ( cyclic(1) ) neighbors_i(:) = mod(neighbors_i(:)+length(1)-grid_shape(1,1), &
01167                                                   length(1))+grid_shape(1,1)
01168             if ( cyclic(2) ) neighbors_j(:) = mod(neighbors_j(:)+length(2)-grid_shape(1,2), &
01169                                                   length(2))+grid_shape(1,2)
01170 
01171             stack(1,stack_entries+1:stack_entries+4) = neighbors_i(:)
01172             stack(2,stack_entries+1:stack_entries+4) = neighbors_j(:)
01173 
01174             stack_entries = stack_entries + 4
01175 
01176          end subroutine add_neigh_to_stack
01177 
01178 !-------------------------------------------------------------------------------------
01179 
01180          subroutine add_boundary_to_stack ( northern, stack_entries )
01181 
01182             use psmile_grid_reduced_gauss
01183 
01184             implicit none
01185 
01186             logical, intent(in)    :: northern
01187             integer, intent(inout) :: stack_entries
01188 
01189             integer :: i
01190 
01191             ! adjust size of stack if necessary
01192             if (size(stack,2) < stack_entries + length(1)) then
01193                stack => psmile_realloc(stack, (/2,size(stack) + &
01194                                                     max(stacksize_inc, length(1))/))
01195             endif
01196 
01197             stack(1,stack_entries+1:stack_entries+length(1)) = &
01198                (/(i, i=grid_shape(1,1), grid_shape(2,1))/)
01199 
01200             if (northern) then
01201                stack(2,stack_entries+1:stack_entries+length(1)) = grid_shape(2,2)
01202             else
01203                stack(2,stack_entries+1:stack_entries+length(1)) = grid_shape(1,2)
01204             endif
01205             
01206             stack_entries = stack_entries + length(1)
01207 
01208          end subroutine add_boundary_to_stack
01209 
01210 !-------------------------------------------------------------------------------------
01211 #ifdef META_POST
01212 
01213          subroutine mp_fill ( s1, s2, s3, s4, overlaps )
01214             use PSMILe, only : point_dble
01215 
01216             Type(point_dble), intent(in) :: s1, s2, s3, s4
01217             Logical, intent(in)          :: overlaps
01218 
01219             Type(point_dble)  :: p1, p2, p3, p4
01220             Character(len=32) :: line_ending
01221 
01222             if (overlaps) then
01223                line_ending = '--cycle withcolor .6*white;'
01224             else
01225                line_ending = '--cycle withcolor .9*white;'
01226             endif
01227             p1 = s1
01228             p2 = s2
01229             p3 = s3
01230             p4 = s4
01231 
01232             p1%x = s1%x + meta_shift_x
01233             p2%x = s2%x + meta_shift_x
01234             p3%x = s3%x + meta_shift_x
01235             p4%x = s4%x + meta_shift_x
01236 
01237             p1%y = s1%y + meta_shift_y
01238             p2%y = s2%y + meta_shift_y
01239             p3%y = s3%y + meta_shift_y
01240             p4%y = s4%y + meta_shift_y
01241 
01242             write(99,'(a6,4(f9.3,a2,f9.3,a5),a27)') 'fill (',p1%x,'u,',p1%y,'u)--(' &
01243                                                             ,p2%x,'u,',p2%y,'u)--(' &
01244                                                             ,p3%x,'u,',p3%y,'u)--(' &
01245                                                             ,p4%x,'u,',p4%y,'u)',trim(line_ending)
01246          end subroutine
01247 
01248 !-------------------------------------------------------------------------------------
01249 
01250          subroutine mp_draw ( s1, s2, is_target )
01251             use PSMILe, only : point_dble
01252             Type(point_dble), intent(in) :: s1, s2
01253             Logical, intent(in)          :: is_target
01254 
01255             Type(point_dble)  :: p1, p2
01256             Character(len=32) :: line_ending
01257 
01258             if (is_target) then
01259                line_ending = '--cycle withcolor black;'
01260             else
01261                line_ending = '--cycle withcolor red;'
01262             endif
01263             p1 = s1
01264             p2 = s2
01265             p1%x = s1%x + meta_shift_x
01266             p2%x = s2%x + meta_shift_x
01267 
01268             p1%y = s1%y + meta_shift_y
01269             p2%y = s2%y + meta_shift_y
01270 
01271             write(99,'(a6,2(f9.3,a2,f9.3,a5),a25)') 'draw (',p1%x,'u,',p1%y,'u)--(' &
01272                                                             ,p2%x,'u,',p2%y,'u)',trim(line_ending)
01273 
01274          end subroutine
01275 #endif
01276       end subroutine psmile_neigh_cells_3d_dble