psmile_mg_next_level_3d_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE: PSMILe_MG_Next_level_3d_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_mg_next_level_3d_dble (            &
00012                            grid_id, lev, nlev,             &
00013                            chmin1, chmin2, chmin3,         &
00014                            chmax1, chmax2, chmax3,         &
00015                            midp1,  midp2,  midp3,          &
00016                            levdim,                         &
00017                            found, loc, range,              &
00018                            coords1, coords2, coords3,      &
00019                            shape, control,                 &
00020                            ijkinc, ijkcoa, ierror)
00021 !
00022 ! !USES:
00023 !
00024       use PRISM_constants
00025 !
00026       use PSMILe, dummy_interface => PSMILe_mg_next_level_3d_dble
00027 
00028       implicit none
00029 !
00030 ! !INPUT PARAMETERS:
00031 !
00032       Integer, Intent (In)            :: grid_id
00033 
00034 !     Info on the component in which the donor
00035 !     should be searched.
00036 
00037       Integer, Intent (In)            :: lev
00038 
00039 !     Level number of the next grid
00040 
00041       Integer, Intent (In)            :: nlev
00042 
00043 !     Total number of levels
00044 
00045       Integer, Intent (In)            :: levdim (ndim_3d)
00046 
00047 !     Dimension of chmin, chmax and midp
00048 
00049       Integer, Intent (In)            :: range (2, ndim_3d)
00050 
00051 !     Dimension of loc and found
00052 
00053       Integer, Intent (In)            :: shape (2, ndim_3d)
00054 
00055 !     Dimension of coordinate arrays
00056 
00057       Double Precision, Intent (In)   :: chmin1 (0:levdim(1), 
00058                                                  0:levdim(2), 
00059                                                  0:levdim(3))
00060       Double Precision, Intent (In)   :: chmin2 (0:levdim(1), 
00061                                                  0:levdim(2), 
00062                                                  0:levdim(3))
00063       Double Precision, Intent (In)   :: chmin3 (0:levdim(1), 
00064                                                  0:levdim(2), 
00065                                                  0:levdim(3))
00066 
00067 !     Minimum of grid coordinates per grid cell
00068 
00069       Double Precision, Intent (In)   :: chmax1 (0:levdim(1), 
00070                                                  0:levdim(2), 
00071                                                  0:levdim(3))
00072       Double Precision, Intent (In)   :: chmax2 (0:levdim(1), 
00073                                                  0:levdim(2), 
00074                                                  0:levdim(3))
00075       Double Precision, Intent (In)   :: chmax3 (0:levdim(1), 
00076                                                  0:levdim(2), 
00077                                                  0:levdim(3))
00078 
00079 !     Maximum of grid coordinates per grid cell
00080 
00081       Double Precision, Intent (In)   ::  midp1 (0:levdim(1), 
00082                                                  0:levdim(2), 
00083                                                  0:levdim(3))
00084       Double Precision, Intent (In)   ::  midp2 (0:levdim(1), 
00085                                                  0:levdim(2), 
00086                                                  0:levdim(3))
00087       Double Precision, Intent (In)   ::  midp3 (0:levdim(1), 
00088                                                  0:levdim(2), 
00089                                                  0:levdim(3))
00090 
00091 !     Mid point of the cell
00092 
00093       Integer, Intent (InOut)         :: found (range(1,1):range(2,1), 
00094                                                 range(1,2):range(2,2), 
00095                                                 range(1,3):range(2,3))
00096 
00097 !     Finest level number on which a grid cell was found for point i,j,k.
00098 !     Level number < -1: Point was not found and
00099 !                        and last level number was (-found(i,j,k))
00100 !     Level number = (-nlev+1): Never found (input value)
00101 
00102       Integer, Intent (InOut)         :: loc   (ndim_3d,               
00103                                                 range(1,1):range(2,1), 
00104                                                 range(1,2):range(2,2), 
00105                                                 range(1,3):range(2,3))
00106 
00107 !     Indices of the grid cell in which the point was found.
00108 
00109       Double Precision, Intent (In)   :: coords1 (shape(1,1):shape(2,1), 
00110                                                   shape(1,2):shape(2,2), 
00111                                                   shape(1,3):shape(2,3))
00112       Double Precision, Intent (In)   :: coords2 (shape(1,1):shape(2,1), 
00113                                                   shape(1,2):shape(2,2), 
00114                                                   shape(1,3):shape(2,3))
00115       Double Precision, Intent (In)   :: coords3 (shape(1,1):shape(2,1), 
00116                                                   shape(1,2):shape(2,2), 
00117                                                   shape(1,3):shape(2,3))
00118 
00119 !     Coordinates to be searched
00120 
00121       Integer, Intent (InOut)         :: control (2, ndim_3d)
00122 
00123 !     Input  value: index range in "coords" to be searched
00124 !     Output value: index range in "coords" for which a cell was found
00125 
00126       Integer, Intent (In)            :: ijkinc (ndim_3d)
00127 
00128 !     Increment
00129 
00130       Integer, Intent (In)            :: ijkcoa (ndim_3d)
00131 
00132 !     Coarsening
00133 !
00134 ! !OUTPUT PARAMETERS:
00135 !
00136       integer, Intent (Out)           :: ierror
00137 
00138 !     Returns the error code of PSMILe_MG_coars_level;
00139 !             ierror = 0 : No error
00140 !             ierror > 0 : Severe error
00141 !
00142 ! !DEFINED PARAMETERS:
00143 !
00144       Integer, parameter              :: ndtry = 64
00145       Integer, parameter              :: ntry1 = 8
00146       Integer, parameter              :: ntry3 = 27
00147       Integer, parameter              :: ntry2 = ndtry - ntry1
00148       Double Precision, parameter     :: remax = 1.0e20
00149 !
00150 ! !LOCAL VARIABLES
00151 !
00152 !     ... for levels
00153 !
00154       Integer                         :: levc
00155 !
00156 !     ... for locations searched
00157 !
00158       Integer                         :: i, j, k
00159       Integer                         :: ibegl, n, nprev
00160       Integer                         :: ifound, nmin(1)
00161       Integer                         :: ijkf(ndim_3d), ijkold(ndim_3d), newijk(ndim_3d)
00162       Integer                         :: ijkctl (ndim_3d, ndtry), ijkdst (ndim_3d, ndtry)
00163       Integer                         :: ijkctl3 (ndim_3d, ndtry)
00164       Double Precision                :: dist (ndtry)
00165       Double Precision                :: xyz (ndim_3d)
00166 !
00167 #ifdef DEBUG
00168 !     ... for locations searched
00169 !
00170       Integer                      :: nfnd (0:4)
00171 #endif /* DEBUG */
00172 !
00173 ! !DESCRIPTION:
00174 !
00175 ! Subroutine "PSMILe_mg_next_level_3d_dble" searches the donor cells
00176 ! on the next MG grid for the subgrid coords by the sending process.
00177 !
00178 !
00179 ! !REVISION HISTORY:
00180 !
00181 !   Date      Programmer   Description
00182 ! ----------  ----------   -----------
00183 ! 03.07.21    H. Ritzdorf  created
00184 !
00185 !EOP
00186 !----------------------------------------------------------------------
00187 !
00188 !  $Id: psmile_mg_next_level_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $
00189 !  $Author: valcke $
00190 !
00191    Character(len=len_cvs_string), save :: mycvs = 
00192        '$Id: psmile_mg_next_level_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
00193 !
00194 !----------------------------------------------------------------------
00195 !
00196 !  Relative control indices for sub-cells of coarse grid cells
00197 !
00198       data ((ijkctl (i, n), i=1,ndim_3d), n = 1,ntry1) &
00199                              / 0, 0, 0,   1, 0, 0,    &
00200                                0, 1, 0,   1, 1, 0,    &
00201                                0, 0, 1,   1, 0, 1,    &
00202                                0, 1, 1,   1, 1, 1 /
00203 !
00204 !  Relative control indices for neighbour cells
00205 !
00206       data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+1,  ntry1+16) &
00207                   /-1,-1,-1,   0,-1,-1,   1,-1,-1,   2,-1,-1,    &
00208                    -1, 0,-1,   0, 0,-1,   1, 0,-1,   2, 0,-1,    &
00209                    -1, 1,-1,   0, 1,-1,   1, 1,-1,   2, 1,-1,    &
00210                    -1, 2,-1,   0, 2,-1,   1, 2,-1,   2, 2,-1/
00211       data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+17, ntry1+28) &
00212                   /-1,-1, 0,   0,-1, 0,   1,-1, 0,   2,-1, 0,    &
00213                    -1, 0, 0,                         2, 0, 0,    &
00214                    -1, 1, 0,                         2, 1, 0,    &
00215                    -1, 2, 0,   0, 2, 0,   1, 2, 0,   2, 2, 0/
00216       data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+29, ntry1+40) &
00217                   /-1,-1, 1,   0,-1, 1,   1,-1, 1,   2,-1, 1,    &
00218                    -1, 0, 1,                         2, 0, 1,    &
00219                    -1, 1, 1,                         2, 1, 1,    &
00220                    -1, 2, 1,   0, 2, 1,   1, 2, 1,   2, 2, 1/
00221       data ((ijkctl (i, n), i=1,ndim_3d), n = ntry1+41, ntry1+56) &
00222                   /-1,-1, 2,   0,-1, 2,   1,-1, 2,   2,-1, 2,    &
00223                    -1, 0, 2,   0, 0, 2,   1, 0, 2,   2, 0, 2,    &
00224                    -1, 1, 2,   0, 1, 2,   1, 1, 2,   2, 1, 2,    &
00225                    -1, 2, 2,   0, 2, 2,   1, 2, 2,   2, 2, 2/
00226 !
00227 !  Relative control indices for computed cells
00228 !
00229       data ((ijkctl3 (i, n), i=1,ndim_3d), n = 1,ntry3) &
00230                  / -1,-1,-1,   0,-1,-1,   1,-1,-1,     &
00231                    -1, 0,-1,   0, 0,-1,   1, 0,-1,     &
00232                    -1, 1,-1,   0, 1,-1,   1, 1,-1,     &
00233                    -1,-1, 0,   0,-1, 0,   1,-1, 0,     &
00234                    -1, 0, 0,   0, 0, 0,   1, 0, 0,     &
00235                    -1, 1, 0,   0, 1, 0,   1, 1, 0,     &
00236                    -1,-1, 1,   0,-1, 1,   1,-1, 1,     &
00237                    -1, 0, 1,   0, 0, 1,   1, 0, 1,     &
00238                    -1, 1, 1,   0, 1, 1,   1, 1, 1 /
00239 !
00240 !----------------------------------------------------------------------
00241 !
00242 !  Initialization
00243 !
00244 #ifdef VERBOSE
00245       print 9990, trim(ch_id), lev, control
00246 
00247       call psmile_flushstd
00248 #endif /* VERBOSE */
00249 !
00250 #ifdef PRISM_ASSERTION
00251 #endif
00252 !
00253       ierror = 0
00254       levc = lev + 1
00255 !
00256 #ifdef DEBUG
00257       nfnd (0) = 0
00258       nfnd (1) = 0
00259       nfnd (2) = 0
00260       nfnd (3) = 0
00261       nfnd (4) = 0
00262 #endif /* DEBUG */
00263 
00264          do k = control(1,3), control(2,3), ijkinc (3)
00265            do j = control(1,2), control (2,2), ijkinc (2)
00266 !
00267 !-----------------------------------------------------------------------
00268 !     Look in the fine parts of the coarse grid cell loc(i,j,k)
00269 !
00270 !     TODO: Index list instead of while loop
00271 !
00272 !     TO TEST:
00273 !         (i)   Compute minimal distance without checking chmin and chmax
00274 !         (ii)  Check chmin and chmax for cell found
00275 !         (iii) If outside, check all cells and mark cells by remax and
00276 !               compute minimum again
00277 !-----------------------------------------------------------------------
00278 !
00279       nprev = 0
00280 
00281       ibegl = control (1, 1)
00282 !
00283 !===> ... Look for the next cell to be controlled
00284 !
00285       do while (ibegl <= control(2,1))
00286 !cdir vector
00287            do i = ibegl, control(2,1), ijkinc(1)
00288            if (found (i,j,k) == levc) go to 11
00289            end do
00290 !
00291         if (i .lt. control(2,1)+ijkinc(1)) then
00292            i = control(2,1)
00293            if (found (i,j,k) == levc) go to 11
00294         endif
00295 !
00296         exit
00297 !
00298 !===> ... Look in the fine parts of cell loc(i,j,k)
00299 !
00300 11         continue
00301 !
00302            ijkf(:) = min (loc(:,i,j,k) * ijkcoa(:), levdim(:))
00303 !
00304 !cdir vector
00305               do n = 1, ntry1
00306               ijkdst (:, n) = max(0, min(ijkf(:) + ijkctl(:,n), levdim(:)))
00307               end do
00308 
00309 !cdir vector
00310               do n = 1, ntry1
00311               if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
00312                                            ijkdst(2,n),         &
00313                                            ijkdst(3,n)) .and.   &
00314                   coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
00315                                            ijkdst(2,n),         &
00316                                            ijkdst(3,n)) .and.   &
00317                   coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
00318                                            ijkdst(2,n),         &
00319                                            ijkdst(3,n)) .and.   &
00320                   coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
00321                                            ijkdst(2,n),         &
00322                                            ijkdst(3,n)) .and.   &
00323                   coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
00324                                            ijkdst(2,n),         &
00325                                            ijkdst(3,n)) .and.   &
00326                   coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
00327                                            ijkdst(2,n),         &
00328                                            ijkdst(3,n))) then
00329                  dist (n) = (coords1(i,j,k) -                                   &
00330                                midp1(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00331                           + (coords2(i,j,k) -                                   &
00332                                midp2(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00333                           + (coords3(i,j,k) -                                   &
00334                                midp3(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2
00335 
00336               else
00337 
00338                  dist (n) = remax
00339 
00340               endif
00341               end do
00342 !
00343 !===> ... Look for the minimum distance
00344 !
00345               nmin = MINLOC (dist(1:ntry1))
00346 #ifdef MINLOCFIX
00347               if (nmin(1).eq.0) nmin=1
00348 #endif /* MINLOCFIX */
00349 !
00350               if (dist(nmin(1)) .ne. remax) then
00351                  found (i,j,k)  = lev
00352                  loc (:,i,j,k) = ijkdst (:,nmin(1))
00353 !
00354 #ifdef DEBUG
00355                  nfnd (1) = nfnd (1) + 1
00356 #endif /* DEBUG */
00357                  go to 95
00358 !
00359               else
00360 !
00361            if (levc == nlev) then
00362               found (i,j,k) = - found(i,j,k)
00363 #ifdef DEBUG
00364               nfnd (0) = nfnd (0) + 1
00365 #endif /* DEBUG */
00366               go to 95
00367            endif
00368         endif
00369 !
00370 !-----------------------------------------------------------------------
00371 !          Look in the neigbourhood of the coarse grid cell
00372 !          loc(i,j,k)
00373 !-----------------------------------------------------------------------
00374 !
00375 !cdir vector
00376               do n = ntry1+1, ndtry
00377               ijkdst (:, n) = max (0, min(ijkf(:) + ijkctl(:,n), levdim(:)))
00378               end do
00379 !
00380               do n = ntry1+1, ndtry
00381               if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
00382                                            ijkdst(2,n),         &
00383                                            ijkdst(3,n)) .and.   &
00384                   coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
00385                                            ijkdst(2,n),         &
00386                                            ijkdst(3,n)) .and.   &
00387                   coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
00388                                            ijkdst(2,n),         &
00389                                            ijkdst(3,n)) .and.   &
00390                   coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
00391                                            ijkdst(2,n),         &
00392                                            ijkdst(3,n)) .and.   &
00393                   coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
00394                                            ijkdst(2,n),         &
00395                                            ijkdst(3,n)) .and.   &
00396                   coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
00397                                            ijkdst(2,n),         &
00398                                            ijkdst(3,n))) then
00399                  dist (n) = (coords1(i,j,k) -                                   &
00400                                midp1(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00401                           + (coords2(i,j,k) -                                   &
00402                                midp2(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00403                           + (coords3(i,j,k) -                                   &
00404                                midp3(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2
00405 
00406               else
00407 
00408                  dist (n) = remax
00409 
00410               endif
00411               end do
00412 !
00413 !===> ... Look for the minimum distance
00414 !
00415            nmin = MINLOC (dist(ntry1+1:ndtry)) + ntry1
00416 #ifdef MINLOCFIX
00417            if (nmin(1).eq.ntry1) nmin=ntry1+1
00418 #endif /* MINLOCFIX */
00419 !
00420            if (dist(nmin(1)) .ne. remax) then
00421               found (i,j,k)  = lev
00422               loc (:,i,j,k) = ijkdst (:,nmin(1))
00423 !
00424 #ifdef DEBUG
00425               nfnd (2) = nfnd (2) + 1
00426 #endif /* DEBUG */
00427 !
00428            else
00429 !
00430               nprev = nprev + 1
00431            endif
00432 !
00433 !===> ... Start index of the next search
00434 !
00435 95         if (i == control(2,1)) exit
00436            ibegl = min (control(2,1), i+ijkinc(1))
00437 
00438          end do ! while
00439 !
00440 !===> ... Have points to be controlled o a coarser level ?
00441 !
00442          if (nprev == 0) go to 20
00443 !
00444 !-----------------------------------------------------------------------
00445 !          Some points were not found.
00446 !          Try to interpolate the fine point a from neighbouring points
00447 !-----------------------------------------------------------------------
00448 !
00449         ibegl = control (1, 1) + ijkinc(1)
00450 !
00451         do while (ibegl <= control(2,1))
00452 !cdir vector
00453               do i = ibegl, control(2,1)-ijkinc(1), ijkinc(1)
00454               if (found (i,j,k)           == levc .and. &
00455                   found (i-ijkinc(1),j,k) == lev  .and. &
00456                   found (i+ijkinc(1),j,k) == lev) go to 211
00457               end do
00458 !
00459            if (i .lt. control(2,1)) then
00460               i = control(2,1) - ijkinc(1)
00461               if (found (i,j,k)           == levc .and. &
00462                   found (i-ijkinc(1),j,k) == lev  .and. &
00463                   found (i+ijkinc(1),j,k) == lev) go to 211
00464            endif
00465 !
00466            exit
00467 !
00468 !===> ... Interpolate cell IJKF and
00469 !         look in the neighbourhood of cell IJKF
00470 !
00471 211        continue
00472 
00473            ijkf (:) = min((loc (:, i-ijkinc(1),j,k) + loc (:, i+ijkinc(1),j,k))/2, &
00474                           levdim(:))
00475 !
00476            ijkold (:) = min (loc(:, i,j,k) * ijkcoa(:), levdim(:))
00477 !
00478            if (ijkold(1) == ijkf(1) .and. &
00479                ijkold(2) == ijkf(2) .and. &
00480                ijkold(3) == ijkf(3)) go to 260
00481 
00482 !cdir vector
00483               do n = 1, ntry3
00484               ijkdst (:, n) = max (0, min(ijkf(:) + ijkctl3(:,n), levdim(:)))
00485               end do
00486 !
00487 !===> ... Compute distances to the cell midpoints
00488 !
00489 !cdir vector
00490               do n = 1, ntry3
00491               if (coords1(i,j,k) >= chmin1(ijkdst(1,n),         &
00492                                            ijkdst(2,n),         &
00493                                            ijkdst(3,n)) .and.   &
00494                   coords2(i,j,k) >= chmin2(ijkdst(1,n),         &
00495                                            ijkdst(2,n),         &
00496                                            ijkdst(3,n)) .and.   &
00497                   coords3(i,j,k) >= chmin3(ijkdst(1,n),         &
00498                                            ijkdst(2,n),         &
00499                                            ijkdst(3,n)) .and.   &
00500                   coords1(i,j,k) <= chmax1(ijkdst(1,n),         &
00501                                            ijkdst(2,n),         &
00502                                            ijkdst(3,n)) .and.   &
00503                   coords2(i,j,k) <= chmax2(ijkdst(1,n),         &
00504                                            ijkdst(2,n),         &
00505                                            ijkdst(3,n)) .and.   &
00506                   coords3(i,j,k) <= chmax3(ijkdst(1,n),         &
00507                                            ijkdst(2,n),         &
00508                                            ijkdst(3,n))) then
00509                  dist (n) = (coords1(i,j,k) -                                   &
00510                                midp1(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00511                           + (coords2(i,j,k) -                                   &
00512                                midp2(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2 &
00513                           + (coords3(i,j,k) -                                   &
00514                                midp3(ijkdst(1,n), ijkdst(2,n), ijkdst(3,n)))**2
00515 
00516               else
00517 
00518                  dist (n) = remax
00519 
00520               endif
00521               end do
00522 !
00523 !===> ... Look for the minimum distance
00524 !
00525            nmin = MINLOC (dist(1:ntry3))
00526 #ifdef MINLOCFIX
00527            if (nmin(1).eq.0) nmin=1
00528 #endif /* MINLOCFIX */
00529 !
00530            if (dist(nmin(1)) .ne. remax) then
00531               found (i,j,k)  = lev
00532               loc (:, i,j,k) = ijkdst (:, nmin(1))
00533 !
00534               nprev = nprev - 1
00535 !
00536 #ifdef DEBUG
00537               nfnd (3) = nfnd (3) + 1
00538 #endif /* DEBUG */
00539         endif
00540 !
00541 !===> ... Start index of the next search
00542 !
00543 260        if (i >= control(2,1)-ijkinc(1)) exit
00544            ibegl = min (control(2,1)-ijkinc(1), i+ijkinc(1))
00545 !
00546          end do ! while
00547 !
00548          if (nprev == 0) go to 20
00549 
00550 !-----------------------------------------------------------------------
00551 !     Some points were still not found.
00552 !     Look again on the next coarser level whether the points
00553 !     are located in other cells.
00554 !-----------------------------------------------------------------------
00555 !
00556       ibegl = control (1, 1)
00557 
00558       do while (ibegl <= control(2,1))
00559 !
00560 !===> ... Look for the next cell to be controlled
00561 !
00562 !cdir vector
00563             do i = ibegl, control(2,1), ijkinc(1)
00564             if (found (i,j,k) == levc) go to 311
00565             end do
00566 !
00567          if (i .lt. control(2,1)+ijkinc(1)) then
00568             i = control(2,1)
00569             if (found (i,j,k) == levc) go to 311
00570          endif
00571 !
00572          exit
00573 !
00574 !===> ... Look in the coarser parts of cell loc(i,j,k)
00575 !
00576 311      continue
00577 !
00578          if (lev+1 < nlev) then
00579             xyz (1) = coords1 (i,j,k)
00580             xyz (2) = coords2 (i,j,k)
00581             xyz (3) = coords3 (i,j,k)
00582 !
00583             call psmile_mg_prev_levels_3d_dble (grid_id, lev, nlev, &
00584                               loc(1,i,j,k), xyz, ifound, newijk)
00585 !
00586             if (ifound .gt. 0) then
00587                loc (:, i,j,k) = newijk (:)
00588 !
00589                found (i,j,k) = lev
00590 #ifdef DEBUG
00591                nfnd (4) = nfnd (4) + 1
00592 #endif /* DEBUG */
00593             else
00594                found (i,j,k) = - found (i,j,k)
00595 #ifdef DEBUG
00596                nfnd (0) = nfnd (0) + 1
00597 #endif /* DEBUG */
00598             endif
00599          else
00600 !           This is the second coarsest level.
00601 !           The coarsest level contains only 1 cell.
00602 !           All points which are not found are outside.
00603 
00604             found (i,j,k) = - found (i,j,k)
00605 #ifdef DEBUG
00606             nfnd (0) = nfnd (0) + 1
00607 #endif /* DEBUG */
00608 
00609          endif
00610 !
00611 !===> ...
00612 !
00613          nprev = nprev - 1
00614          if (nprev == 0) go to 20
00615 !
00616 !===> ... Start index of the next search
00617 !
00618          if (i == control(2,1)) exit
00619 
00620          ibegl = min (control(2,1), i+ijkinc(1))
00621       end do ! while
00622 !
00623 20    end do ! j
00624       end do ! k
00625 !
00626 !===> All done
00627 !
00628 #ifdef DEBUG
00629       print 9970, trim(ch_id), lev, nfnd, ijkinc
00630 9970  format (1x, a, ': PSMILe_mg_next_level_3d_dble: lev =', i3, &
00631               ': fnd =', 5i5, ', ijkinc ', 3i4)
00632 #endif /* DEBUG */
00633 
00634 #ifdef VERBOSE
00635       print 9980, trim(ch_id), lev
00636 
00637       call psmile_flushstd
00638 #endif /* VERBOSE */
00639 !
00640 !  Formats:
00641 !
00642 9990 format (1x, a, ': PSMILe_mg_next_level_3d_dble: level', i3, &
00643                     ', control', 6i6)
00644 9980 format (1x, a, ': PSMILe_mg_next_level_3d_dble: eof, level', i3)
00645 
00646       end subroutine PSMILe_mg_next_level_3d_dble