psmile_mg_cells_2d_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 ! !ROUTINE: PSMILe_MG_cells_2d_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_mg_cells_2d_dble ( grid_id,       &
00012                            search_grid_type,              &
00013                            found, loc, loc_fnd_shape,     &
00014                            tgt_corners_x, tgt_corners_y,  &
00015                            tgt_corner_shape, control,     &
00016                            grid_valid_shape, ipart,       &
00017                            src_corner_shape, nbr_corners, &
00018                            src_corners_x, src_corners_y,  &
00019                            chmin1, chmax1,                &
00020                            chmin2, chmax2,                &
00021                            tol, ierror)
00022 !
00023 ! !USES:
00024 !
00025       use PRISM_constants
00026       use psmile_grid, only : common_grid_range
00027       use PSMILe, dummy_interface => PSMILe_mg_cells_2d_dble
00028 #ifdef DEBUG_TRACE
00029       use psmile_debug_trace
00030 #endif
00031 
00032       implicit none
00033 !
00034 ! !INPUT PARAMETERS:
00035 !
00036       Integer, Intent (In)            :: grid_id
00037 
00038       Integer, Intent (In)            :: search_grid_type
00039 
00040       Integer, Intent (In)            :: loc_fnd_shape (2, ndim_3d)
00041 
00042 !     Dimension of loc and found
00043 
00044       Integer, Intent (In)            :: tgt_corner_shape (2, ndim_3d)
00045 
00046 !     Dimension of target corner arrays
00047 
00048       Integer, Intent (InOut)         :: found (loc_fnd_shape(1,1):loc_fnd_shape(2,1), 
00049                                                 loc_fnd_shape(1,2):loc_fnd_shape(2,2), 
00050                                                 loc_fnd_shape(1,3):loc_fnd_shape(2,3))
00051 
00052 !     Finest level number on which a grid cell was found for point i,j,k.
00053 !     Input :
00054 !     Level number < -1: Point was not found and
00055 !                        and last level number was (-found(i,j,k))
00056 !     Level number = -(nlev+1): Never found
00057 !     Output :
00058 !     found(i,j,k) = +1: Point (i,j,k) is located on a point
00059 !
00060 !     found(i,j,k) = -1: Point (i,j,k) is located in a cell
00061 !
00062 !     Otherwise, not contained in the cell grid.
00063 
00064       Integer, Intent (InOut)         :: loc   (ndim_2d,               
00065                                                 loc_fnd_shape(1,1):loc_fnd_shape(2,1), 
00066                                                 loc_fnd_shape(1,2):loc_fnd_shape(2,2), 
00067                                                 loc_fnd_shape(1,3):loc_fnd_shape(2,3))
00068 
00069 !     Indices of the grid cell in which the point was found.
00070 !     The indices are relative to "tgt_corner_shape".
00071 
00072       Double Precision, Intent (In)   :: tgt_corners_x (                               
00073                                           tgt_corner_shape(1,1):tgt_corner_shape(2,1), 
00074                                           tgt_corner_shape(1,2):tgt_corner_shape(2,2), 
00075                                           tgt_corner_shape(1,3):tgt_corner_shape(2,3))
00076       Double Precision, Intent (In)   :: tgt_corners_y (                               
00077                                           tgt_corner_shape(1,1):tgt_corner_shape(2,1), 
00078                                           tgt_corner_shape(1,2):tgt_corner_shape(2,2), 
00079                                           tgt_corner_shape(1,3):tgt_corner_shape(2,3))
00080 
00081 !     Target corners to be searched
00082 
00083       Integer, Intent (In)            :: control (2, ndim_3d)
00084 
00085 !     Input  value: index range in "tgt_corners" to be searched
00086 !     Output value: index range in "tgt_corners" for which a cell was found
00087 
00088       Integer,          Intent (In)   :: grid_valid_shape (2,2)
00089 
00090 !     Specifies the valid block shape for the "1d"-dimensional block
00091 !     Dimension of chmin and chmax
00092 
00093       Integer,          Intent (In)   :: ipart
00094 
00095       Integer,          Intent (In)   :: src_corner_shape (2,2)
00096 
00097       Integer,          Intent (In)   :: nbr_corners
00098 
00099       Double Precision, Intent (In)   :: src_corners_x (                                
00100                                           src_corner_shape (1,1):src_corner_shape(2,1), 
00101                                           src_corner_shape (1,2):src_corner_shape(2,2), 
00102                                           nbr_corners)
00103 
00104       Double Precision, Intent (In)   :: src_corners_y (                                
00105                                           src_corner_shape (1,1):src_corner_shape(2,1), 
00106                                           src_corner_shape (1,2):src_corner_shape(2,2), 
00107                                           nbr_corners)
00108 
00109       Double Precision, Intent (InOut)  :: chmin1 (grid_valid_shape(1,1):    
00110                                                    grid_valid_shape(2,1)+1 , 
00111                                                    grid_valid_shape(1,2):    
00112                                                    grid_valid_shape(2,2)+1)
00113       Double Precision, Intent (InOut)  :: chmin2 (grid_valid_shape(1,1):    
00114                                                    grid_valid_shape(2,1)+1,  
00115                                                    grid_valid_shape(1,2):    
00116                                                    grid_valid_shape(2,2)+1)
00117 !
00118 !     Minimum value of bounding box of coordinate.
00119 ! chmin (grid_valid_shape(1)+i) contains the minimal value of bounding box of
00120 !                   cell (src_corners_x(grid_valid_shape(1)+i)  :
00121 !                         src_corners_x(grid_valid_shape(1)+i+1))
00122 !
00123 ! The first box
00124 ! chmin (grid_valid_shape(1)-1) contains the minimal value of
00125 !                                        the bounding box from
00126 !                    the virtual cell (src_corners_x(grid_valid_shape(1)):
00127 !                                               bounding box of
00128 !                                               corner control volume)
00129 ! The last box
00130 ! chmin (grid_valid_shape(2)) contains the minimal value of
00131 !                                      the bounding box from
00132 !                    the virtual cell (src_corners_x(grid_valid_shape(2)):
00133 !                                               bounding box of
00134 !                                               corner control volume)
00135 !
00136       Double Precision, Intent (InOut)  :: chmax1 (grid_valid_shape(1,1):    
00137                                                    grid_valid_shape(2,1)+1 , 
00138                                                    grid_valid_shape(1,2):    
00139                                                    grid_valid_shape(2,2)+1)
00140       Double Precision, Intent (InOut)  :: chmax2 (grid_valid_shape(1,1):    
00141                                                    grid_valid_shape(2,1)+1,  
00142                                                    grid_valid_shape(1,2):    
00143                                                    grid_valid_shape(2,2)+1)
00144 !
00145 !     Maximum value of bounding box of coordinate.
00146 !
00147       Double Precision, Intent (In)   :: tol
00148 !
00149 !     Absolute tolerance for search of "identical" points
00150 !     TOL >= 0.0
00151 !
00152 ! !OUTPUT PARAMETERS:
00153 !
00154       integer, Intent (Out)           :: ierror
00155 
00156 !     Returns the error code of PSMILE_mg_cells_2d_dble;
00157 !             ierror = 0 : No error
00158 !             ierror > 0 : Severe error
00159 !
00160 ! !DEFINED PARAMETERS:
00161 !
00162 !  NREFD = Number of method coordinates to be controlled.
00163 !  val_direct  = Code for locations which should to be directly transferred
00164 !                to the destination process.
00165 !  val_coupler = Code for locations which should to be transferred
00166 !                to the coupler process.
00167 
00168       Integer, Parameter              :: nrefd = 3
00169 !
00170       Integer, Parameter              :: val_direct  =  1
00171       Integer, Parameter              :: val_coupler = -1
00172 
00173       Integer, Parameter              :: temp_len = 256
00174 !
00175 !     ... for levels
00176 !
00177       Integer, Parameter              :: lev = 1
00178 !
00179 ! !LOCAL VARIABLES
00180 !
00181 !
00182 !     Total number of levels
00183 !
00184       Integer                         :: nlev
00185 !
00186 !     Period of a cyclic coordinate.
00187 !
00188       Double Precision                :: period (ndim_2d)
00189 !
00190 !     ... for locations searched
00191 !
00192       Type (Corner_Block), Pointer    :: corner_pointer
00193       Integer                         :: control_ (2, ndim_3d) !this is actually just a copy of
00194                                                                !control which is slightly modified
00195       Integer                         :: leni, lenij
00196       Integer                         :: i, j, k, l
00197       Integer                         :: ibegl, iendl
00198       Integer                         :: ic(temp_len, ndim_2d)
00199       Integer                         :: ijk(3)
00200 !
00201 !     ... For locations controlled
00202 !
00203       Integer                         :: irange(temp_len, 2, ndim_2d)
00204       Integer                         :: ioff ! Offset in Gaussred data structure
00205 !
00206 !     ... for search
00207 !
00208       Integer                         :: ii, jj, nc
00209       Integer                         :: list (temp_len*2)
00210 !   
00211 !     ... For cyclic cells
00212 !
00213       Double Precision                :: shifted (ndim_2d)
00214 !
00215 #ifdef DEBUG
00216 !     ... for locations searched
00217 !
00218       Integer                         :: nfnd (0:3)
00219 #endif /* DEBUG */
00220 !
00221 ! !DESCRIPTION:
00222 !
00223 ! Subroutine "PSMILe_mg_cells_2d_dble" searches the corresponding cells
00224 ! of the grid for the subgrid corners of the sending process.
00225 !
00226 ! TODO: Improve !!
00227 !       Use tol(erance) for overlap tests. Same tol has to be used later in
00228 !       PSMILe_Neigh_cells_*
00229 !
00230 ! !REVISION HISTORY:
00231 !
00232 !   Date      Programmer   Description
00233 ! ----------  ----------   -----------
00234 ! 06.08.29    R. Redler    created
00235 !
00236 !EOP
00237 !----------------------------------------------------------------------
00238 !
00239 !  $Id: psmile_mg_cells_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $
00240 !  $Author: hanke $
00241 !
00242    Character(len=len_cvs_string), save :: mycvs = 
00243        '$Id: psmile_mg_cells_2d_dble.F90 2927 2011-01-28 14:04:12Z hanke $'
00244 !----------------------------------------------------------------------
00245 !
00246 !  Initialization
00247 !
00248       ierror = 0
00249       corner_pointer => Grids(grid_id)%corner_pointer
00250 
00251       nlev   = Grids(grid_id)%nlev
00252       period = common_grid_range(2,1:2) - common_grid_range(1,1:2)
00253 
00254 #ifdef VERBOSE
00255       print 9990, trim(ch_id), lev, control
00256 
00257       call psmile_flushstd
00258 #endif /* VERBOSE */
00259 !
00260 #ifdef DEBUGX
00261       print *, ' tgt_corner_shape ', tgt_corner_shape
00262       print *, ' loc_fnd_shape ', loc_fnd_shape
00263 
00264       do k = loc_fnd_shape(1,3), loc_fnd_shape(2,3)
00265          do j = loc_fnd_shape(1,2), loc_fnd_shape(2,2)
00266             do i = loc_fnd_shape(1,1), loc_fnd_shape(2,1)
00267                print *, ' loc and fnd ', i, j, k, ':', &
00268                                          loc(1,i,j,k), &
00269                                          loc(2,i,j,k), &
00270                                          found(i,j,k)
00271             enddo
00272          enddo
00273       enddo
00274 #endif
00275 !
00276 #ifdef PRISM_ASSERTION
00277       if (tol < 0.0) then
00278          call psmile_assert (__FILE__, __LINE__, &
00279                              "tol must be >= 0.0")
00280       endif
00281 #endif
00282 !
00283 ! create temporary memory (control (2,1)-control(1,1)+1)*ndtry
00284 ! and do n =1,nref before i-loop ?!?
00285 !
00286 ! ? instead of controlling all nref, determine location in cell once
00287 !
00288 ! ? one-dimensional Loop instead of i,j loop
00289 !
00290 #ifdef DEBUG
00291       nfnd (:) = 0
00292 #endif /* DEBUG */
00293 
00294 #ifdef DEBUG_TRACE
00295       if (loc_fnd_shape(1,1) <= ictl_ind(1) .and. loc_fnd_shape(2,1) >= ictl_ind(1) .and. &
00296           loc_fnd_shape(1,2) <= ictl_ind(2) .and. loc_fnd_shape(2,2) >= ictl_ind(2)) then
00297           print *, 'a) ipart treated is ', ipart
00298           print *, 'corners', tgt_corners_x(ictl_ind(1),ictl_ind(2),control(1,3)), &
00299                               tgt_corners_y(ictl_ind(1),ictl_ind(2),control(1,3))
00300           print 8990, ictl_ind (1:ndim_2d), &
00301                       found ( ictl_ind(1),ictl_ind(2),control(1,3)), &
00302                       loc (:, ictl_ind(1),ictl_ind(2),control(1,3))
00303       endif
00304 #endif /* DEBUG_TRACE */
00305 
00306 #ifdef DEBUGX
00307       print *, 'begin ---'
00308       do k = control (1, 3), control (2,3)   
00309          do j = control(1,2), control (2,2) 
00310             do i = control(1,1), control (2,1) 
00311                print *, 'i, j, k, tgt_corners_x, tgt_corners_y, locx, locy, found', i,j,k, &
00312                     tgt_corners_x (i,j,k),  tgt_corners_y (i,j,k), &
00313                     loc(1,i,j,k), loc(2,i,j,k), found (i,j,k)
00314             end do
00315          end do
00316       end do
00317 #endif
00318 
00319       control_ = control
00320       ! in case the target grid is not of type PRISM_Gaussreduced_regvrt
00321       ! the upper bound of the first two dimension is reduced by one, because
00322       ! we are only interested in the "lower left" corners of the cells, because
00323       ! they are our reference points
00324       if ( search_grid_type /= PRISM_Gaussreduced_regvrt ) &
00325          control_(2,1:2) = control (2,1:2) - 1
00326 
00327       if ( Grids(grid_id)%pole_covered ) then
00328          leni  =   src_corner_shape(2,1) - src_corner_shape(1,1) + 1
00329          lenij = ( src_corner_shape(2,2) - src_corner_shape(1,2) + 1 ) * leni
00330       endif
00331 !
00332 !===> Define chmin and chmax
00333 !
00334       if ( ipart == 1 ) then
00335          chmin1 =  999.0
00336          chmax1 = -999.0
00337          chmin2 =  999.0
00338          chmax2 = -999.0
00339 
00340          do j = grid_valid_shape(1,2), grid_valid_shape(2,2)
00341             do i = grid_valid_shape(1,1), grid_valid_shape(2,1)
00342                chmin1(i,j) = minval(src_corners_x(i,j,:))
00343                chmax1(i,j) = maxval(src_corners_x(i,j,:))
00344                chmin2(i,j) = minval(src_corners_y(i,j,:))
00345                chmax2(i,j) = maxval(src_corners_y(i,j,:))
00346             enddo
00347          enddo
00348          !
00349          ! Extend to pole if necessary
00350          !
00351          if ( Grids(grid_id)%pole_covered ) then
00352 
00353             do l = 1, size(corner_pointer%pole_array)
00354 
00355                ijk = psmile_transform_index_1d_to_3d( &
00356                   corner_pointer%pole_array(l), &
00357                   corner_pointer%corner_shape)
00358 
00359                if ( chmax2(ijk(1),ijk(2)) > 0.0 ) then
00360                   chmax2(ijk(1),ijk(2)) =  90.0
00361                else if ( chmin2(ijk(1),ijk(2)) < 0.0 ) then
00362                   chmin2(ijk(1),ijk(2)) = -90.0
00363                endif
00364 
00365             enddo
00366 
00367          endif
00368 
00369       endif
00370 !
00371 !===> Determine location
00372 !
00373       do k = control_ (1, 3), control_ (2,3)
00374          do j = control_(1,2), control_ (2,2)
00375 !
00376 !-----------------------------------------------------------------------
00377 !        Look in the fine parts of the coarse grid cell loc(i)
00378 !-----------------------------------------------------------------------
00379 !
00380             ibegl = control_ (1, 1)
00381 !
00382 !===> ... Look for the next cell to be controlled
00383 !
00384 loop_i:     do while (ibegl <= control_(2,1))
00385 
00386                nc = 0
00387                i = ibegl
00388                do while (nc < temp_len .and. i <= control_(2,1))
00389                   ibegl = i
00390 !cdir vector
00391                   do i = ibegl, min (ibegl+(temp_len*2-1)-nc, control_(2,1))
00392                      if (found(i, j, k) == lev) then
00393                         nc = nc + 1
00394                         list (nc) = i
00395                      endif
00396                   end do
00397                end do
00398 
00399                if (nc == 0) exit loop_i ! exit while loop (ibegl < ...)
00400 
00401                if (nc < temp_len) then
00402                   ibegl = control_(2,1) + 1
00403                else  if (nc > temp_len) then
00404                   nc = temp_len
00405                   ibegl = list (temp_len+1)
00406                else
00407                   ibegl = list (temp_len) + 1
00408                endif
00409 !
00410 !===> ... Look in the real cell (ic(i, 1),ic(i, 2))
00411 !
00412                ic (1:nc, 1) = loc (1, list(1:nc), j, k)
00413                ic (1:nc, 2) = loc (2, list(1:nc), j, k)
00414 
00415 #ifdef PRISM_ASSERTION
00416 !cdir vector
00417                do i = 1, nc
00418                   if (ic(i,1) < grid_valid_shape(1,1) .or. &
00419                       ic(i,1) > grid_valid_shape(2,1) .or. &
00420                       ic(i,2) < grid_valid_shape(1,2) .or. &
00421                       ic(i,2) > grid_valid_shape(2,2) ) then
00422                      exit
00423                   endif
00424                end do
00425 
00426                if (i <= nc) then
00427                   print *, 'i, ic', i, ic(i, :), found(list(i),j,k)
00428                   print *, 'grid_valid_shape', grid_valid_shape
00429                   call psmile_assert (__FILE__, __LINE__, &
00430                        'wrong index')
00431                endif
00432 #endif
00433                do i = 1, nc
00434                   irange (i, 1, :) = max (ic(i,:)-1, grid_valid_shape (1, :))
00435                   irange (i, 2, :) = min (ic(i,:)+1, grid_valid_shape (2, :))
00436                end do
00437 !
00438 !cdir vector
00439                do i = ibegl, control_(2,1)
00440                   if (found(i, j, k) /= lev) exit
00441                end do
00442 !
00443                iendl = i - 1
00444 !
00445 !===> ... Look in the real cell (ic)
00446 !
00447                do i = 1, nc
00448 
00449                   ! move target longitudes into the common local longitude range
00450                   shifted (1) = tgt_corners_x (list(i),j,k)
00451                   do while (tgt_corners_x (list(i),j,k) - period (1) >= chmin1(ic(i,1),ic(i,2)))
00452                      shifted (1) = shifted (1) - period(1)
00453                   enddo
00454                   do while (tgt_corners_x (list(i),j,k) + period (1) < chmax1(ic(i,1),ic(i,2)))
00455                      shifted (1) = shifted (1) + period(1)
00456                   enddo
00457 
00458                   ! we do not allow cyclic corner handling in latitude direction
00459                   shifted (2) = tgt_corners_y (list(i),j,k)
00460 
00461                   do jj = irange (i, 1, 2), irange (i, 2, 2)
00462                      do ii = irange (i, 1, 1), irange (i, 2, 1)
00463 
00464                         if ( shifted(1) >= chmin1(ii,jj) .and.   &
00465                              shifted(2) >= chmin2(ii,jj) .and.   &
00466                              shifted(1) <= chmax1(ii,jj) .and.   &
00467                              shifted(2) <= chmax2(ii,jj)) then
00468 
00469                              ic(i,1) = ii
00470                              ic(i,2) = jj
00471                           
00472                         endif
00473 
00474                      end do ! ii
00475                   end do ! jj
00476 
00477                   ii = ic(i,1)
00478                   jj = ic(i,2)
00479 
00480                   if ( chmin1(ic(i,1),ic(i,2)) > shifted(1) ) ii = ii - 1
00481                   if ( chmin2(ic(i,1),ic(i,2)) > shifted(2) ) jj = jj - 1
00482 
00483                   if ( ii <=  irange (i, 2, 1) .and. ii >= irange (i, 1, 1) .and. &
00484                        jj <=  irange (i, 2, 2) .and. jj >= irange (i, 1, 2) ) then
00485 
00486                      loc (1,list(i),j,k) = ii
00487                      loc (2,list(i),j,k) = jj
00488 
00489 #ifdef DEBUG
00490                      nfnd (1) = 0
00491                      nfnd (2) = nfnd (2) + 1
00492 #endif /* DEBUG */
00493                      found (list(i),j,k) = val_coupler
00494 
00495                   else
00496 !
00497 !              Cell has no overlap with any other cell in the search grid,
00498 !              mark point as 0
00499 #ifdef DEBUG
00500                      nfnd (3) = nfnd (3) + 1
00501 !                    nfnd (2) = nfnd (2) - 1
00502 #endif /* DEBUG */
00503 !Todo                found (i,j,k) = 0
00504 
00505                   endif
00506 
00507                end do ! i
00508 !
00509 !===> ... Start index of the next search
00510 !
00511             ibegl = iendl + 2
00512             end do loop_i ! while
00513 !
00514          end do ! j
00515       end do ! k
00516 !
00517 !===> Check that location found are in the valid range, otherwise mark point as not found.
00518 !     Enforce to send all data through the coupler.
00519 !
00520       do k = loc_fnd_shape (1,3), loc_fnd_shape (2,3)
00521          do j = loc_fnd_shape(1,2), loc_fnd_shape (2,2)
00522             do i = loc_fnd_shape(1,1), loc_fnd_shape (2,1)
00523                if ( abs(found(i,j,k)) == 1 ) found(i,j,k) = val_coupler
00524                if ( loc(1,i,j,k) < grid_valid_shape(1,1) .or. &
00525                     loc(1,i,j,k) > grid_valid_shape(2,1) .or. &
00526                     loc(2,i,j,k) < grid_valid_shape(1,2) .or. &
00527                     loc(2,i,j,k) > grid_valid_shape(2,2) )    &
00528                     found(i,j,k) = -(nlev+1)
00529             end do
00530          end do
00531       end do
00532 !
00533 !===> ... Check whether missed point a located in a pole cell 
00534 !
00535 ! TODO: store nbr_corner_face in data structure in prism_set_corners.
00536 
00537       if ( Grids(grid_id)%pole_covered ) then
00538 
00539          do k = loc_fnd_shape (1,3), loc_fnd_shape (2,3)
00540             do j = loc_fnd_shape(1,2), loc_fnd_shape (2,2)
00541                do i = loc_fnd_shape(1,1), loc_fnd_shape (2,1)
00542 
00543                   if ( abs(found(i,j,k)) /= lev ) then
00544 
00545                      do l = 1, size(corner_pointer%pole_array)
00546 
00547                         ijk = psmile_transform_index_1d_to_3d( &
00548                            corner_pointer%pole_array(l), &
00549                            corner_pointer%corner_shape)
00550 
00551                         if ( tgt_corners_y(i,j,k) >= chmin2(ijk(1),ijk(2)) .and. &
00552                              tgt_corners_y(i,j,k) <= chmax2(ijk(1),ijk(2)) ) then
00553 
00554                            found(i,j,k) = lev
00555                            loc(1,i,j,k) = ijk(1)
00556                            loc(2,i,j,k) = ijk(2)
00557 #ifdef DEBUG
00558                            nfnd(0) = nfnd(0) - 1
00559                            nfnd(2) = nfnd(2) + 1
00560 #endif
00561                         endif
00562 
00563                      enddo ! l-loop
00564 
00565                   endif ! found
00566 
00567                enddo
00568             enddo
00569          enddo
00570 
00571       endif ! pole_coverd
00572 !
00573 !===> Additional check for empty cells. The initial search with mg_method only searched
00574 !     for the minimum corners. Here we check whether one of the other corners falls into
00575 !     the valid area and in this case mark the cell as found.
00576 !
00577       if ( search_grid_type == PRISM_Gaussreduced_regvrt ) then
00578 
00579          ioff = (loc_fnd_shape(2,1)-loc_fnd_shape(1,1)+1)/2
00580          ibegl = loc_fnd_shape(1,1)
00581          iendl = loc_fnd_shape(1,1) + ioff - 1
00582 
00583 #ifdef DEBUG_TRACE
00584       if (loc_fnd_shape(1,1) <= ictl_ind(1) .and. loc_fnd_shape(2,1) >= ictl_ind(1) .and. &
00585           loc_fnd_shape(1,2) <= ictl_ind(2) .and. loc_fnd_shape(2,2) >= ictl_ind(2)) then
00586           print *, ' ipart ', ipart, ' ictl found ', &
00587               found ( ictl_ind(1)     ,ictl_ind(2),control(1,3)), &
00588               found ( ictl_ind(1)+ioff,ictl_ind(2),control(1,3))
00589 
00590            print *, ' ... with corners '
00591 
00592            print *, 'min corners', tgt_corners_x(ictl_ind(1),ictl_ind(2),control(1,3)), &
00593                                    tgt_corners_y(ictl_ind(1),ictl_ind(2),control(1,3))
00594 
00595            print *, 'max corners', tgt_corners_x(ictl_ind(1)+ioff,ictl_ind(2),control(1,3)), &
00596                                    tgt_corners_y(ictl_ind(1)+ioff,ictl_ind(2),control(1,3))
00597        endif
00598 #endif
00599          do k = control_ (1, 3), control_ (2,3)
00600             do j = control_(1,2), control_ (2,2) ! should be "1" for PRISM_Gaussreduced_regvrt
00601                do i = ibegl, iendl
00602                   if ( abs(found(i,j,k)) > 1 ) then
00603                      if ( abs(found(ioff+i,j,k)) == 1 ) then
00604                         loc(1,i,j,k) = loc(1,ioff+i,j,k)
00605                         loc(2,i,j,k) = loc(2,ioff+i,j,k)
00606                         found(i,j,k) = found(ioff+i,j,k)
00607 !rr                     else if ( abs(found(ioff+i-1,j,k)) == 1 ) then
00608 !rr                        loc(1,i,j,k) = loc(1,ioff+i-1,j,k)
00609 !rr                        loc(2,i,j,k) = loc(2,ioff+i-1,j,k)
00610 !rr                        found(i,j,k) = found(ioff+i-1,j,k)
00611 !rr                     else if ( abs(found(i+1,j,k)) == 1 ) then
00612 !rr                        loc(1,i,j,k) = loc(1,i+1,j,k)
00613 !rr                        loc(2,i,j,k) = loc(2,i+1,j,k)
00614 !rr                        found(i,j,k) = found(i+1,j,k)
00615 !                      else
00616 ! 
00617 !                         ! if no corner of the target cell is within a source cell
00618 !                         ! but the corner was in a bounding box of the source data of this process
00619 !                         if (abs(found(i,j,k)) <= nlev) then
00620 ! 
00621 !                            ! it might have happend that the target cell has an overlap with a
00622 !                            ! source cell of this process, but none of its corners is in one
00623 !                            ! of these source cells, in this case we try to find source courners
00624 !                            ! which are with this target cell
00625 ! 
00626 !                            ! get the target cell (for gaussreduced grids the cells are rectangular)
00627 !                            ! I do not really know why this is correct, I just copied it from the DEBUG_TRACE
00628 !                            ! part above
00629 !                            cell(1,1) = tgt_corners_x(i,j,k)
00630 !                            cell(2,1) = tgt_corners_y(i,j,k)
00631 !                            cell(1,2) = tgt_corners_x(i+ioff,j,k)
00632 !                            cell(2,2) = tgt_corners_y(i+ioff,j,k)
00633 ! 
00634 !                            ! check whether there are source corners within the target cells
00635 ! 
00636 !                     outer: do a = grid_valid_shape(1,1), grid_valid_shape(2,1)
00637 !                               do b = grid_valid_shape(1,2), grid_valid_shape(2,2)
00638 !                                  do c = 1, nbr_corners
00639 !                                     if (inside_bb(cell(:,1:2), (/src_corners_x(a,b,c), src_corners_y(a,b,c)/))) then
00640 !                                        loc(1,i,j,k) = a
00641 !                                        loc(2,i,j,k) = b
00642 !                                        found(i,j,k) = -1
00643 !                                        exit outer
00644 !                                     endif
00645 !                                  enddo
00646 !                               enddo
00647 !                            enddo outer
00648 !                         endif
00649                      endif
00650                   endif
00651                enddo
00652             enddo
00653          enddo
00654 
00655          do k = control_ (1, 3), control_ (2,3)
00656             do j = control_(1,2), control_ (2,2)
00657                do i = iendl+1, loc_fnd_shape(2,1)
00658                   found(i,j,k) = -(nlev+1)
00659                enddo
00660             enddo
00661          enddo
00662 
00663       else
00664 
00665          do k = control_ (1, 3), control_ (2,3)
00666             do j = control_(1,2), control_ (2,2)
00667                do i = control_(1,1), control_ (2,1)
00668 
00669                   if ( abs(found(i,j,k)) > 1 ) then
00670                      if ( abs(found(i+1,j,k)) == 1 ) then
00671                         found(i,j,k) = found(i+1,j,k)
00672                         loc(1,i,j,k) = loc(1,i+1,j,k)
00673                         loc(2,i,j,k) = loc(2,i+1,j,k)
00674                      else if ( abs(found(i,j+1,k)) == 1 ) then
00675                         found(i,j,k) = found(i,j+1,k)
00676                         loc(1,i,j,k) = loc(1,i,j+1,k)
00677                         loc(2,i,j,k) = loc(2,i,j+1,k)
00678                      else if ( abs(found(i+1,j+1,k)) == 1 ) then
00679                         found(i,j,k) = found(i+1,j+1,k)
00680                         loc(1,i,j,k) = loc(1,i+1,j+1,k)
00681                         loc(2,i,j,k) = loc(2,i+1,j+1,k)
00682 !                      else if (abs(found(i,j,k)) <= nlev) then
00683 ! 
00684 !                         ! get the target cell
00685 !                         cell(1,1) = tgt_corners_x(i  ,j  ,k)
00686 !                         cell(2,1) = tgt_corners_y(i  ,j  ,k)
00687 !                         cell(1,2) = tgt_corners_x(i+1,j  ,k)
00688 !                         cell(2,2) = tgt_corners_y(i+1,j  ,k)
00689 !                         cell(1,3) = tgt_corners_x(i+1,j+1,k)
00690 !                         cell(2,3) = tgt_corners_y(i+1,j+1,k)
00691 !                         cell(1,4) = tgt_corners_x(i  ,j+1,k)
00692 !                         cell(2,4) = tgt_corners_y(i  ,j+1,k)
00693 !                         ! create a bounding box
00694 !                         bounding_box(1,1) = minval(cell(1,1:4))
00695 !                         bounding_box(2,1) = minval(cell(2,1:4))
00696 !                         bounding_box(1,2) = maxval(cell(1,1:4))
00697 !                         bounding_box(2,2) = maxval(cell(2,1:4))
00698 ! 
00699 !                         ! check whether there are source corners within the target cells
00700 !                         ! and use them instead (see handling of gaussreduced grids for more comments)
00701 !                  outer2: do a = grid_valid_shape(1,1), grid_valid_shape(2,1)
00702 !                            do b = grid_valid_shape(1,2), grid_valid_shape(2,2)
00703 !                               do c = 1, nbr_corners
00704 !                                  if (inside_bb(bounding_box, (/src_corners_x(a,b,c), src_corners_y(a,b,c)/))) then
00705 !                                     ! an additional check for the exact cell is still missing
00706 !                                     loc(1,i,j,k) = a
00707 !                                     loc(2,i,j,k) = b
00708 !                                     found(i,j,k) = -1
00709 !                                     exit outer2
00710 !                                  endif
00711 !                               enddo
00712 !                            enddo
00713 !                         enddo outer2
00714                      endif
00715                   endif
00716                enddo
00717             enddo
00718          enddo
00719 
00720          do k = control_ (1, 3), control_ (2,3)
00721             j = control_ (2,2)+1
00722             do i = control_(1,1), control_ (2,1)+1
00723                found(i,j,k) = -(nlev+1)
00724             enddo
00725             i = control_ (2,1)+1
00726             do j = control_(1,2), control_ (2,2)+1
00727                found(i,j,k) = -(nlev+1)
00728             enddo
00729          enddo
00730 
00731       endif
00732 !
00733 !===> All done
00734 !
00735 #ifdef DEBUG_TRACE
00736       if (loc_fnd_shape(1,1) <= ictl_ind(1) .and. loc_fnd_shape(2,1) >= ictl_ind(1) .and. &
00737           loc_fnd_shape(1,2) <= ictl_ind(2) .and. loc_fnd_shape(2,2) >= ictl_ind(2)) then
00738           print *, 'b) ipart treated is ', ipart
00739           print *, 'tgt_corners', tgt_corners_x(ictl_ind(1),ictl_ind(2),control(1,3)), &
00740                                   tgt_corners_y(ictl_ind(1),ictl_ind(2),control(1,3))
00741           print 8990, ictl_ind (1:ndim_2d), &
00742                       found ( ictl_ind(1),ictl_ind(2),control(1,3)), &
00743                       loc (:, ictl_ind(1),ictl_ind(2),control(1,3))
00744 8990  format (1x, '### psmile_mg_cells_2d_dble: ictl_ind', 2i5, &
00745                   '; found, loc ', i3, 1x, 2i6)
00746       endif
00747 #endif /* DEBUG_TRACE */
00748 
00749 #ifdef DEBUGX
00750       print *, 'end ---'
00751         do k = control (1, 3), control (2,3)
00752             do j = control(1,2), control (2,2)
00753                do i = control(1,1), control (2,1)
00754                print *, 'i, j, k, tgt_corners_x, tgt_corners_y, locx, locy, found', i,j,k, &
00755                         tgt_corners_x (i,j,k),  tgt_corners_y (i,j,k),  &
00756                         loc(1,i,j,k), loc(2,i,j,k), found (i,j,k)
00757                end do
00758             end do
00759         end do
00760 #endif
00761 
00762 #ifdef DEBUG
00763       print 9970, trim(ch_id), lev, nfnd
00764 9970  format (1x, a, ': psmile_mg_cells_2d_dble: lev =', i3, &
00765               ': not :', i5, ', dir: ', i6, ', fnd :', i6, i5)
00766 #endif /* DEBUG */
00767 
00768 #ifdef VERBOSE
00769       print 9980, trim(ch_id), lev
00770 
00771       call psmile_flushstd
00772 #endif /* VERBOSE */
00773 !
00774 !  Formats:
00775 !
00776 9990 format (1x, a, ': psmile_mg_cells_2d_dble: level', i3, &
00777                     ', control', 6i6)
00778 9980 format (1x, a, ': psmile_mg_cells_2d_dble: eof, level', i3)
00779 
00780       contains
00781 
00782          logical function inside_bb(bb, p)
00783             double precision :: bb(2,2), p(2)
00784 
00785             if (p(1) >= bb(1,1) .and. p(1) <= bb(1,2) .and. &
00786                 p(2) >= bb(2,1) .and. p(2) <= bb(2,2)) then
00787 
00788                inside_bb = .true.
00789             else
00790                inside_bb = .false.
00791             endif
00792 
00793          end function inside_bb
00794 
00795       end subroutine psmile_mg_cells_2d_dble
Generated on 18 Mar 2011 for Oasis4 by  doxygen 1.6.1