psmile_bbcells_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_Bbcells_3d_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_bbcells_3d_dble (method_id, array, shape,     &
00012                            range, corners, corner_shape, nbr_corners, &
00013                            chmin,   chmax,  midp,      levdim,        &
00014                            cyclic,  period, ierror)
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019 !
00020       use PSMILe, dummy_interface => PSMILe_Bbcells_3d_dble
00021 
00022       implicit none
00023 !
00024 ! !INPUT PARAMETERS:
00025 !
00026       Integer,           Intent (In)  :: method_id
00027 
00028 !     Method id to access halo information
00029 
00030       Integer,           Intent (In)  :: shape (2, ndim_3d)
00031 
00032 !     Dimension of (method/grid) coordinate "array"
00033 
00034       Double Precision,  Intent (In)  :: array(shape(1,1):shape(2,1), 
00035                                                shape(1,2):shape(2,2), 
00036                                                shape(1,3):shape(2,3))
00037 !     Coordinates of the method
00038 
00039       Integer,           Intent (In)  :: range (2, ndim_3d)
00040 
00041 !     Definition of the subgrid for which the (bounding) cells
00042 !     should be computed.
00043 
00044       Integer,           Intent (In)  :: corner_shape (2, ndim_3d)
00045 
00046 !     Dimension of the array "corners" containing the corner points
00047 !     of the control volume (for a single coordinate).
00048 
00049       Integer,           Intent (In)  :: nbr_corners
00050 
00051 !     Number of corners per control volume
00052 
00053       Double Precision,  Intent (In)  :: corners (corner_shape(1,1): 
00054                                                   corner_shape(2,1), 
00055                                                   corner_shape(1,2): 
00056                                                   corner_shape(2,2), 
00057                                                   corner_shape(1,3): 
00058                                                   corner_shape(2,3), 
00059                                                   nbr_corners)
00060 
00061 !     Definintion of the subgrid for which the (bounding)
00062 !     cells should be computed.
00063 
00064       Integer,           Intent (In)  :: levdim (ndim_3d)
00065 
00066 !     Dimension of output arrays "chmin", "chmax" and "midp".
00067 !     For testing.
00068 !
00069       Logical,           Intent (In)  :: cyclic
00070 !
00071       Double Precision,  Intent (In)  :: period
00072 !
00073 ! !OUTPUT PARAMETERS:
00074 !
00075       Double Precision,  Intent (Out) :: chmin (range(1,1)-1:range(2,1), 
00076                                                 range(1,2)-1:range(2,2), 
00077                                                 range(1,3)-1:range(2,3))
00078 !
00079 !     Minimum value of bounding box of coordinate.
00080 !
00081       Double Precision,  Intent (Out) :: chmax (range(1,1)-1:range(2,1), 
00082                                                 range(1,2)-1:range(2,2), 
00083                                                 range(1,3)-1:range(2,3))
00084 !
00085 !     Maximum value of bounding box of coordinate.
00086 !
00087       Double Precision,  Intent (Out) :: midp  (range(1,1)-1:range(2,1), 
00088                                                 range(1,2)-1:range(2,2), 
00089                                                 range(1,3)-1:range(2,3))
00090 !
00091 !     Mid point of bounding box of coordinate.
00092 !
00093       Integer,           Intent (Out) :: ierror
00094 
00095 !     Returns the error code of PSMILe_Bbcells_3d_dble
00096 !             ierror = 0 : No error
00097 !             ierror > 0 : Severe error
00098 !
00099 ! !DEFINED PARAMETERS:
00100 !
00101 ! n_corners_3d = Number of corners in a 3d block
00102 ! r_nbr        = 1 / Number of corners per method cell.
00103 !
00104       Integer,          Parameter  :: n_corners_3d = 8
00105       Double Precision, Parameter  :: r_nbr = 0.125d0
00106 !
00107 ! !LOCAL VARIABLES
00108 !
00109       Integer                      :: i, ibeg, iend
00110       Integer                      :: j, jbeg, jend
00111       Integer                      :: k, kbeg, kend
00112 !
00113       Integer                      :: iv, jv, kv
00114 !
00115 !  Auxiliary vectors for edges and corners
00116 !
00117       Integer                      :: iin (n_corners_3d), jin (n_corners_3d)
00118       Integer                      :: kin (n_corners_3d)
00119       Integer                      :: ivt (n_corners_3d), jvt (n_corners_3d)
00120       Integer                      :: kvt (n_corners_3d)
00121       Integer                      :: ic  (n_corners_3d), jc  (n_corners_3d)
00122       Integer                      :: kc  (n_corners_3d)
00123       Integer                      :: n
00124 !
00125       Integer, Allocatable         :: lchmin (:, :), lchmax (:, :)
00126 !
00127 !  Auxiliary vectors for dimensions
00128 !
00129       Integer                      :: shape_2d (2, ndim_2d)
00130       Integer                      :: range_2d (2, ndim_2d)
00131       Integer                      :: corner_shape_2d (2, ndim_2d)
00132       Integer                      :: levdim_2d (ndim_2d)
00133 !
00134 ! Error parameters
00135 !
00136       Integer, Parameter           :: nerrp = 2
00137       Integer                      :: ierrp (nerrp)
00138 !
00139 ! !DESCRIPTION:
00140 !
00141 ! Subroutine "PSMILe_Bbcells_3d_dble" computes the bounding box for each cell
00142 ! of the 3-dimensional coordinate defined by "array".
00143 !
00144 ! chmin (i, j, k) returns the minimal value of bounding box of
00145 !              method cell ((i,j,k  ), (i+1,j,k  ), (i,j+1,k  ), (i+1,j+1,k  ),
00146 !                           (i,j,k+1), (i+1,j,k+1), (i,j+1,k+1), (i+1,j+1,k+1))
00147 !
00148 ! The first box
00149 ! chmin (range(1,1)-1,j,k) returns the minimal value of the bounding box
00150 !                        of virtual cell
00151 !       (array(range(1,1), j, k  ), array(range(1,1), j+1, k  ),
00152 !        array(range(1,1), j, k+1), array(range(1,1), j+1, k+1),
00153 !        bounding box of corner control volume of (range(1,1),j,k)
00154 !
00155 ! The last box
00156 ! chmin (range(2,1),  j,k) returns the minimal value of the bounding box
00157 !                        of virtual cell
00158 !       (array(range(2,1), j, k  ), array(range(2,1), j+1, k  ),
00159 !        array(range(2,1), j, k+1), array(range(2,1), j+1, k+1),
00160 !        bounding box of corner control volume of (range(2,1),j,k)
00161 !
00162 ! Corresponding values are returned for other virtual cells at block
00163 ! border.
00164 !
00165 ! TODO: (*) Cyclic uebergeben und in der Berechnung der virtuellen
00166 !           Zellen verwenden.
00167 !       (*) Eigentlich braeuchte man die method coordinates of
00168 !           neighbouring blocks.
00169 !       (*) Sonst braucht man nicht beide Coordinaten,
00170 !           um die Entscheidung lokal zu treffen.
00171 !       (*) Fuer psmile_mg_method_... braucht man eigentlich auch
00172 !           die Corners the Methoden Cellen; dh. ein Output Array
00173 !           cell_corners (range(1,1)-1:range(2,1), &
00174 !                         range(1,2)-1:range(2,2), 8)
00175 !           Das bedeutet aber, dass man die Methoden Coordinaten
00176 !           kopieren muesste.  Das wuerde dann die Berechnung der
00177 !           chmin, chmax vereinfachen.
00178 !           Das muss ich nochmal gruendlich ueberlegen.
00179 !
00180 ! !REVISION HISTORY:
00181 !
00182 !   Date      Programmer   Description
00183 ! ----------  ----------   -----------
00184 ! 03.06.25    H. Ritzdorf  created
00185 !
00186 !EOP
00187 !----------------------------------------------------------------------
00188 !
00189 !  $Id: psmile_bbcells_3d_dble.F90 2687 2010-10-28 15:15:52Z coquart $
00190 !  $Author: coquart $
00191 !
00192    Character(len=len_cvs_string), save :: mycvs = 
00193        '$Id: psmile_bbcells_3d_dble.F90 2687 2010-10-28 15:15:52Z coquart $'
00194 !
00195 !----------------------------------------------------------------------
00196 !
00197 #ifdef VERBOSE
00198       print 9990, trim(ch_id)
00199 
00200       call psmile_flushstd
00201 #endif /* VERBOSE */
00202 !
00203 #ifdef PRISM_ASSERTION
00204       if (range(1,1) < corner_shape(1,1) .or. &
00205           range(2,1) > corner_shape(2,1) .or. &
00206           range(1,2) < corner_shape(1,2) .or. &
00207           range(2,2) > corner_shape(2,2) .or. &
00208           range(1,3) < corner_shape(1,3) .or. &
00209           range(2,3) > corner_shape(2,3)) then
00210          print *, trim(ch_id), "PSMILe_bbcells_3d_dble: range", range
00211          print *, trim(ch_id), "PSMILe_bbcells_3d_dble: corner_shape", &
00212                   corner_shape
00213 
00214          call psmile_assert (__FILE__, __LINE__, &
00215               "range must be in corner_shape")
00216       endif
00217 !
00218       if (minval(levdim(:)-(range(2,:)-range(1,:)+2)) < 0) then
00219          print *, trim(ch_id), "PSMILe_bbcells_3d_dble: levdim", levdim
00220          print *, trim(ch_id), "PSMILe_bbcells_3d_dble: range ", range
00221 
00222          call psmile_assert (__FILE__, __LINE__, &
00223               "levdim < range(2,:)-range(1,:)+2")
00224       endif
00225 #endif
00226 !
00227 !  Initialization
00228 !
00229       ierror = 0
00230 !
00231       ibeg = range (1,1)
00232       iend = range (2,1)
00233 !
00234       jbeg = range (1,2)
00235       jend = range (2,2)
00236 !
00237       kbeg = range (1,3)
00238       kend = range (2,3)
00239 !
00240 !=======================================================================
00241 !     Special case: singular array
00242 !=======================================================================
00243 !
00244       if (range(1,3) == range(2,3)) then
00245 !
00246         call psmile_bbcells_2d_dble (                &
00247             array (:, :, kbeg),                         &
00248             shape (:, 1:ndim_2d), range (:, 1:ndim_2d), &
00249             corner_shape (:, 1:ndim_2d),                &
00250             chmin (:, :, kbeg),                         &
00251             chmax (:, :, kbeg),                         &
00252             midp  (:, :, kbeg),                         &
00253             levdim(1:ndim_2d), period, ierror)
00254         if (ierror > 0) return
00255 !
00256 !      ... Copy bounding box to virtual face (default)
00257 !
00258          chmin (:, jbeg:jend, kbeg-1) = chmin (:, jbeg:jend, kbeg)
00259          chmax (:, jbeg:jend, kbeg-1) = chmax (:, jbeg:jend, kbeg)
00260          midp  (:, jbeg:jend, kbeg-1) = midp  (:, jbeg:jend, kbeg)
00261 !
00262 !         ... Allocate work vector for corner volumes
00263 !
00264          Allocate (lchmin (ibeg:iend-1, jbeg:jend-1), &
00265                    lchmax (ibeg:iend-1, jbeg:jend-1), STAT = ierror)
00266          if (ierror > 0) then
00267             ierrp (1) = ierror
00268             ierrp (2) = (iend-ibeg) * (jend-jbeg) * 2
00269 
00270             ierror = PRISM_Error_Alloc
00271 
00272             call psmile_error ( ierror, 'lchmin, lchmax', &
00273                                 ierrp, 2, __FILE__, __LINE__ )
00274             return
00275          endif
00276 !
00277             do j = jbeg, jend-1
00278 !cdir vector
00279                do i = ibeg, iend-1
00280                lchmin (i,j) = min (minval (corners(i,   j,   kbeg, :)), &
00281                                    minval (corners(i+1, j,   kbeg, :)), &
00282                                    minval (corners(i,   j+1, kbeg, :)), &
00283                                    minval (corners(i+1, j+1, kbeg, :)))
00284                enddo
00285             enddo
00286 !
00287             do j = jbeg, jend-1
00288 !cdir vector
00289                do i = ibeg, iend-1
00290                lchmax (i,j) = max (maxval (corners(i,   j,   kbeg, :)), &
00291                                    maxval (corners(i+1, j,   kbeg, :)), &
00292                                    maxval (corners(i,   j+1, kbeg, :)), &
00293                                    maxval (corners(i+1, j+1, kbeg, :)))
00294                enddo
00295             end do
00296 !
00297 !          ... Control and Create final boundig box
00298 !              The decision of the first and last box on line
00299 !              is made using next interior box.
00300 !
00301 ! jv  = Index of       virtual method cell - in j-direction
00302 ! jc  = Index of corresponding corner cell - in j-direction
00303 ! j   = Index of last interior method cell - in j-direction
00304 !
00305          j  = jbeg
00306          jv = jbeg - 1
00307 !        jc = jbeg = j
00308 !
00309 !cdir vector
00310             do i = ibeg, iend-1
00311             if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00312                 chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00313                chmin (i,jv,kbeg)   = array (i,j,kbeg)
00314                chmax (i,jv,kbeg-1) = array (i,j,kbeg)
00315 !
00316                chmin (i,jv,kbeg-1) = minval (corners(i,j,kbeg, :))
00317                chmax (i,jv,kbeg  ) = maxval (corners(i,j,kbeg, :))
00318             endif
00319             enddo
00320 !
00321          j  = jend - 1
00322          jv = jend
00323 !        jc = jend = jv
00324 !
00325 !cdir vector
00326             do i = ibeg, iend-1
00327             if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00328                 chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00329                chmin (i,jv,kbeg)   = array (i,jv,kbeg)
00330                chmax (i,jv,kbeg-1) = array (i,jv,kbeg)
00331 !
00332                chmin (i,jv,kbeg-1) = minval (corners(i,jv,kbeg, :))
00333                chmax (i,jv,kbeg  ) = maxval (corners(i,jv,kbeg, :))
00334             endif
00335             enddo
00336 !
00337          i  = ibeg
00338          iv = ibeg - 1
00339 !        ic = ibeg = i
00340 !
00341 !cdir vector
00342             do j = jbeg, jend-1
00343             if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00344                 chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00345                chmin (iv,j,kbeg)   = array (i,j,kbeg)
00346                chmax (iv,j,kbeg-1) = array (i,j,kbeg)
00347 !
00348                chmin (iv,j,kbeg-1) = minval (corners(i,j,kbeg, :))
00349                chmax (iv,j,kbeg  ) = maxval (corners(i,j,kbeg, :))
00350             endif
00351             enddo
00352 !
00353          i  = iend - 1
00354          iv = iend
00355 !        ic = iend = iv
00356 !
00357 !cdir vector
00358             do j = jbeg, jend-1
00359             if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00360                 chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00361                chmin (iv,j,kbeg)   = array (iv,j,kbeg)
00362                chmax (iv,j,kbeg-1) = array (iv,j,kbeg)
00363 !
00364                chmin (iv,j,kbeg-1) = minval (corners(iv,j,kbeg, :))
00365                chmax (iv,j,kbeg  ) = maxval (corners(iv,j,kbeg, :))
00366             endif
00367             enddo
00368 !
00369 !        Corners:
00370 !
00371          i = ibeg; iv = ibeg-1
00372          j = jbeg; jv = jbeg-1
00373 !
00374          if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00375              chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00376             chmin (iv,jv,kbeg)   = array (ibeg,jbeg,kbeg)
00377             chmax (iv,jv,kbeg-1) = array (ibeg,jbeg,kbeg)
00378 !
00379             chmin (iv,jv,kbeg-1) = minval (corners(ibeg,jbeg,kbeg, :))
00380             chmax (iv,jv,kbeg  ) = maxval (corners(ibeg,jbeg,kbeg, :))
00381          endif
00382 !
00383          i = iend-1; iv = iend
00384          j = jbeg;   jv = jbeg-1
00385 !
00386          if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00387              chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00388             chmin (iv,jv,kbeg)   = array (iend,jbeg,kbeg)
00389             chmax (iv,jv,kbeg-1) = array (iend,jbeg,kbeg)
00390 !
00391             chmin (iv,jv,kbeg-1) = minval (corners(iend,jbeg,kbeg, :))
00392             chmax (iv,jv,kbeg  ) = maxval (corners(iend,jbeg,kbeg, :))
00393          endif
00394 !
00395          i = ibeg;   iv = ibeg-1
00396          j = jend-1; jv = jend
00397 !
00398          if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00399              chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00400             chmin (iv,jv,kbeg)   = array (ibeg,jend,kbeg)
00401             chmax (iv,jv,kbeg-1) = array (ibeg,jend,kbeg)
00402 !
00403             chmin (iv,jv,kbeg-1) = minval (corners(ibeg,jend,kbeg, :))
00404             chmax (iv,jv,kbeg  ) = maxval (corners(ibeg,jend,kbeg, :))
00405          endif
00406 !
00407          i = iend-1; iv = iend
00408          j = jend-1; jv = jend
00409 !
00410          if (chmin (i,j,kbeg) >= minval (corners(i,j,kbeg, :)) .and. &
00411              chmax (i,j,kbeg) <= maxval (corners(i,j,kbeg, :))) then
00412             chmin (iv,jv,kbeg)   = array (iend,jend,kbeg)
00413             chmax (iv,jv,kbeg-1) = array (iend,jend,kbeg)
00414 !
00415             chmin (iv,jv,kbeg-1) = minval (corners(iend,jend,kbeg, :))
00416             chmax (iv,jv,kbeg  ) = maxval (corners(iend,jend,kbeg, :))
00417          endif
00418 !
00419 !        Lower and upper Face
00420 !
00421             do j = jbeg, jend-1
00422 !cdir vector
00423                do i = ibeg, iend-1
00424                if (chmin (i,j, kbeg) >= minval (corners(i,j, kbeg, :)) .and. &
00425                    chmax (i,j, kbeg) <= maxval (corners(i,j, kbeg, :))) then
00426                   chmin (i,j, kbeg)   = (chmin(i,j, kbeg) + chmax(i,j, kbeg)) * 0.5
00427                   chmax (i,j, kbeg-1) = chmin (i,j, kbeg)
00428 !
00429                   chmin (i,j, kbeg-1) = lchmin (i,j)
00430                   chmax (i,j, kbeg  ) = lchmax (i,j)
00431                endif
00432                end do
00433             end do
00434 !
00435          Deallocate (lchmin, lchmax)
00436 !
00437 !----------------------------------------------------------------------
00438 !    Singular in j-direction
00439 !----------------------------------------------------------------------
00440 !
00441       else if (range (1,2) == range (2,2)) then
00442 !
00443         shape_2d (:, 1) = shape (:, 1)
00444         shape_2d (:, 2) = shape (:, 3)
00445 !
00446         range_2d (:, 1) = range (:, 1)
00447         range_2d (:, 2) = range (:, 3)
00448 !
00449         corner_shape_2d (:, 1) = corner_shape (:, 1)
00450         corner_shape_2d (:, 2) = corner_shape (:, 3)
00451 !
00452         levdim_2d (1) = levdim (1)
00453         levdim_2d (2) = levdim (3)
00454 !
00455         call psmile_bbcells_2d_dble (                &
00456             array (:, :, kbeg),                         &
00457             shape (:, 1:ndim_2d), range (:, 1:ndim_2d), &
00458             corner_shape (:, 1:ndim_2d),                &
00459             chmin (:, :, kbeg),                         &
00460             chmax (:, :, kbeg),                         &
00461             midp  (:, :, kbeg),                         &
00462             levdim(1:ndim_2d), period, ierror)
00463         if (ierror > 0) return
00464 !
00465 !      ... Copy bounding box to virtual face (default)
00466 !
00467          chmin (:, jbeg-1, kbeg:kend) = chmin (:, jbeg, kbeg:kend)
00468          chmax (:, jbeg-1, kbeg:kend) = chmax (:, jbeg, kbeg:kend)
00469          midp  (:, jbeg-1, kbeg:kend) = midp  (:, jbeg, kbeg:kend)
00470 !
00471 !         ... Allocate work vector for corner volumes
00472 !
00473          Allocate (lchmin (ibeg:iend-1, kbeg:kend-1), &
00474                    lchmax (ibeg:iend-1, kbeg:kend-1), STAT = ierror)
00475          if (ierror > 0) then
00476             ierrp (1) = ierror
00477             ierrp (2) = (iend-ibeg) * (kend-kbeg) * 2
00478 
00479             ierror = PRISM_Error_Alloc
00480 
00481             call psmile_error ( ierror, 'lchmin, lchmax', &
00482                                 ierrp, 2, __FILE__, __LINE__ )
00483             return
00484          endif
00485 !
00486             do k = kbeg, kend-1
00487 !cdir vector
00488                do i = ibeg, iend-1
00489                lchmin (i,k) = min (minval (corners(i,   jbeg, k  , :)), &
00490                                    minval (corners(i+1, jbeg, k  , :)), &
00491                                    minval (corners(i,   jbeg, k+1, :)), &
00492                                    minval (corners(i+1, jbeg, k+1, :)))
00493                enddo
00494             enddo
00495 !
00496             do k = kbeg, kend-1
00497 !cdir vector
00498                do i = ibeg, iend-1
00499                lchmax (i,k) = max (maxval (corners(i,   jbeg, k  , :)), &
00500                                    maxval (corners(i+1, jbeg, k  , :)), &
00501                                    maxval (corners(i,   jbeg, k+1, :)), &
00502                                    maxval (corners(i+1, jbeg, k+1, :)))
00503                enddo
00504             end do
00505 !
00506 !          ... Control and Create final boundig box
00507 !              The decision of the first and last box on line
00508 !              is made using next interior box.
00509 !
00510 ! kv  = Index of       virtual method cell - in k-direction
00511 ! kc  = Index of corresponding corner cell - in k-direction
00512 ! k   = Index of last interior method cell - in k-direction
00513 !
00514          k  = kbeg
00515          kv = kbeg - 1
00516 !        kc = kbeg = k
00517 !
00518 !cdir vector
00519             do i = ibeg, iend-1
00520             if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00521                 chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00522                chmin (i,jbeg,  kv) = array (i,jbeg,k)
00523                chmax (i,jbeg-1,kv) = array (i,jbeg,k)
00524 !
00525                chmin (i,jbeg-1,kv) = minval (corners(i,jbeg,k, :))
00526                chmax (i,jbeg,  kv) = maxval (corners(i,jbeg,k, :))
00527             endif
00528             enddo
00529 !
00530          k  = kend - 1
00531          kv = kend
00532 !        kc = kend = kv
00533 !
00534 !cdir vector
00535             do i = ibeg, iend-1
00536             if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00537                 chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00538                chmin (i,jbeg,  kv) = array (i,jbeg,kv)
00539                chmax (i,jbeg-1,kv) = array (i,jbeg,kv)
00540 !
00541                chmin (i,jbeg-1,kv) = minval (corners(i,jbeg,kv, :))
00542                chmax (i,jbeg,  kv) = maxval (corners(i,jbeg,kv, :))
00543             endif
00544             enddo
00545 !
00546          i  = ibeg
00547          iv = ibeg - 1
00548 !        ic = ibeg = i
00549 !
00550 !cdir vector
00551             do k = kbeg, kend-1
00552             if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00553                 chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00554                chmin (iv,jbeg  ,k) = array (i,jbeg,k)
00555                chmax (iv,jbeg-1,k) = array (i,jbeg,k)
00556 !
00557                chmin (iv,jbeg-1,k) = minval (corners(i,jbeg,k, :))
00558                chmax (iv,jbeg,  k) = maxval (corners(i,jbeg,k, :))
00559             endif
00560             enddo
00561 !
00562          i  = iend - 1
00563          iv = iend
00564 !        ic = iend = iv
00565 !
00566 !cdir vector
00567             do k = kbeg, kend-1
00568             if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00569                 chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00570                chmin (iv,jbeg,  k) = array (iv,jbeg,k)
00571                chmax (iv,jbeg-1,k) = array (iv,jbeg,k)
00572 !
00573                chmin (iv,jbeg-1,k) = minval (corners(iv,jbeg,k, :))
00574                chmax (iv,jbeg,  k) = maxval (corners(iv,jbeg,k, :))
00575             endif
00576             enddo
00577 !
00578 !        Corners
00579 !
00580          i = ibeg; iv = ibeg-1
00581          k = kbeg; kv = kbeg-1 ! ; kc = kbeg
00582 !
00583          if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00584              chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00585             chmin (iv,jbeg  ,kv) = array (ibeg,jbeg,kbeg)
00586             chmax (iv,jbeg-1,kv) = array (ibeg,jbeg,kbeg)
00587 !
00588             chmin (iv,jbeg-1,kv) = minval (corners(ibeg,jbeg,kbeg, :))
00589             chmax (iv,jbeg  ,kv) = maxval (corners(ibeg,jbeg,kbeg, :))
00590          endif
00591 !
00592          i = iend-1; iv = iend
00593          k = kbeg;   kv = kbeg-1 ! ; kc = kbeg
00594 !
00595          if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00596              chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00597             chmin (iv,jbeg,  kv) = array (iend,jbeg,kbeg)
00598             chmax (iv,jbeg-1,kv) = array (iend,jbeg,kbeg)
00599 !
00600             chmin (iv,jbeg-1,kv) = minval (corners(iend,jbeg,kbeg, :))
00601             chmax (iv,jbeg  ,kv) = maxval (corners(iend,jbeg,kbeg, :))
00602          endif
00603 !
00604          i = ibeg;   iv = ibeg-1
00605          k = kend-1; kv = kend   ! ; kc = kend
00606 !
00607          if (chmin (i,jend,k) >= minval (corners(i,jend,k, :)) .and. &
00608              chmax (i,jend,k) <= maxval (corners(i,jend,k, :))) then
00609             chmin (iv,jend,  kv) = array (ibeg,jend,kend)
00610             chmax (iv,jend-1,kv) = array (ibeg,jend,kend)
00611 !
00612             chmin (iv,jend-1,kv) = minval (corners(ibeg,jend,kend, :))
00613             chmax (iv,jend  ,kv) = maxval (corners(ibeg,jend,kend, :))
00614          endif
00615 !
00616          i = iend-1; iv = iend
00617          k = kend-1; kv = kend ! ; kc = kend
00618 !
00619          if (chmin (i,jend,k) >= minval (corners(i,jend,k, :)) .and. &
00620              chmax (i,jend,k) <= maxval (corners(i,jend,k, :))) then
00621             chmin (iv,jend  ,kv) = array (iend,jend,kend)
00622             chmax (iv,jend-1,kv) = array (iend,jend,kend)
00623 !
00624             chmin (iv,jend-1,kv) = minval (corners(iend,jend,kend, :))
00625             chmax (iv,jend  ,kv) = maxval (corners(iend,jend,kend, :))
00626          endif
00627 !
00628 !        Lower and upper Face
00629 !
00630             do k = kbeg, kend-1
00631 !cdir vector
00632                do i = ibeg, iend-1
00633                if (chmin (i,jbeg,k) >= minval (corners(i,jbeg,k, :)) .and. &
00634                    chmax (i,jbeg,k) <= maxval (corners(i,jbeg,k, :))) then
00635                   chmin (i,jbeg,   k) = (chmin(i,jbeg,k) + chmax(i,jbeg,k)) * 0.5
00636                   chmax (i,jbeg-1, k) = chmin (i,jbeg,k)
00637 !
00638                   chmin (i,jbeg-1, k) = lchmin (i,k)
00639                   chmax (i,jbeg,   k) = lchmax (i,k)
00640                endif
00641                end do
00642             end do
00643 !
00644          Deallocate (lchmin, lchmax)
00645 !
00646 !----------------------------------------------------------------------
00647 !    Singular in i-direction
00648 !----------------------------------------------------------------------
00649 !
00650       else if (range (1,1) == range (2,1)) then
00651 !
00652         call psmile_bbcells_2d_dble (                &
00653             array (:, :, kbeg),                         &
00654             shape (:, 1:ndim_2d), range (:, 1:ndim_2d), &
00655             corner_shape (:, 1:ndim_2d),                &
00656             chmin (:, :, kbeg),                         &
00657             chmax (:, :, kbeg),                         &
00658             midp  (:, :, kbeg),                         &
00659             levdim(1:ndim_2d), period, ierror)
00660         if (ierror > 0) return
00661 !
00662 !      ... Copy bounding box to virtual face (default)
00663 !
00664          chmin (ibeg-1, jbeg:jend, kbeg:kend) = chmin (ibeg, jbeg:jend, kbeg:kend)
00665          chmax (ibeg-1, jbeg:jend, kbeg:kend) = chmax (ibeg, jbeg:jend, kbeg:kend)
00666          midp  (ibeg-1, jbeg:jend, kbeg:kend) = midp  (ibeg, jbeg:jend, kbeg:kend)
00667 !
00668 !         ... Allocate work vector for corner volumes
00669 !
00670          Allocate (lchmin (jbeg:jend-1, kbeg:kend-1), &
00671                    lchmax (jbeg:jend-1, kbeg:kend-1), STAT = ierror)
00672          if (ierror > 0) then
00673             ierrp (1) = ierror
00674             ierrp (2) = (jend-jbeg) * (kend-kbeg) * 2
00675 
00676             ierror = PRISM_Error_Alloc
00677 
00678             call psmile_error ( ierror, 'lchmin, lchmax', &
00679                                 ierrp, 2, __FILE__, __LINE__ )
00680             return
00681          endif
00682 !
00683             do k = kbeg, kend-1
00684 !cdir vector
00685                do j = jbeg, jend-1
00686                lchmin (j,k) = min (minval (corners(ibeg, j  , k  , :)), &
00687                                    minval (corners(ibeg, j+1, k  , :)), &
00688                                    minval (corners(ibeg, j  , k+1, :)), &
00689                                    minval (corners(ibeg, j+1, k+1, :)))
00690                enddo
00691             enddo
00692 !
00693             do k = kbeg, kend-1
00694 !cdir vector
00695                do j = jbeg, jend-1
00696                lchmax (j,k) = max (maxval (corners(ibeg, j  , k  , :)), &
00697                                    maxval (corners(ibeg, j+1, k  , :)), &
00698                                    maxval (corners(ibeg, j  , k+1, :)), &
00699                                    maxval (corners(ibeg, j+1, k+1, :)))
00700                enddo
00701             end do
00702 !
00703 !          ... Control and Create final boundig box
00704 !              The decision of the first and last box on line
00705 !              is made using next interior box.
00706 !
00707 ! kv  = Index of       virtual method cell - in k-direction
00708 ! kc  = Index of corresponding corner cell - in k-direction
00709 ! k   = Index of last interior method cell - in k-direction
00710 !
00711          k  = kbeg
00712          kv = kbeg - 1
00713 !        kc = kbeg = k
00714 !
00715 !cdir vector
00716             do j = jbeg, jend-1
00717             if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00718                 chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00719                chmin (ibeg,  j,kv) = array (ibeg,j,k)
00720                chmax (ibeg-1,j,kv) = array (ibeg,j,k)
00721 !
00722                chmin (ibeg-1,j,kv) = minval (corners(ibeg,j,k, :))
00723                chmax (ibeg,  j,kv) = maxval (corners(ibeg,j,k, :))
00724             endif
00725             enddo
00726 !
00727          k  = kend - 1
00728          kv = kend
00729 !        kc = kend = kv
00730 !
00731 !cdir vector
00732             do j = jbeg, jend-1
00733             if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00734                 chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00735                chmin (ibeg  ,j,kv) = array (ibeg,j,kv)
00736                chmax (ibeg-1,j,kv) = array (ibeg,j,kv)
00737 !
00738                chmin (ibeg-1,j,kv) = minval (corners(ibeg,j,kv, :))
00739                chmax (ibeg,  j,kv) = maxval (corners(ibeg,j,kv, :))
00740             endif
00741             enddo
00742 !
00743          j  = jbeg
00744          jv = jbeg - 1
00745 !        jc = jbeg = j
00746 !
00747 !cdir vector
00748             do k = kbeg, kend-1
00749             if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00750                 chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00751                chmin (ibeg  ,jv,k) = array (ibeg,j,k)
00752                chmax (ibeg-1,jv,k) = array (ibeg,j,k)
00753 !
00754                chmin (ibeg-1,jv,k) = minval (corners(ibeg,j,k, :))
00755                chmax (ibeg  ,jv,k) = maxval (corners(ibeg,j,k, :))
00756             endif
00757             enddo
00758 !
00759          j  = jend - 1
00760          jv = jend
00761 !        jc = jend = jv
00762 !
00763 !cdir vector
00764             do k = kbeg, kend-1
00765             if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00766                 chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00767                chmin (ibeg,  jv,k) = array (ibeg,jv,k)
00768                chmax (ibeg-1,jv,k) = array (ibeg,jv,k)
00769 !
00770                chmin (ibeg-1,jv,k) = minval (corners(ibeg,jv,k, :))
00771                chmax (ibeg,  jv,k) = maxval (corners(ibeg,jv,k, :))
00772             endif
00773             enddo
00774 !
00775 !        Corners
00776 !
00777          j = jbeg; jv = jbeg-1
00778          k = kbeg; kv = kbeg-1 ! ; kc = kbeg
00779 !
00780          if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00781              chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00782             chmin (ibeg,  jv,kv) = array (ibeg,jbeg,kbeg)
00783             chmax (ibeg-1,jv,kv) = array (ibeg,jbeg,kbeg)
00784 !
00785             chmin (ibeg-1,jv,kv) = minval (corners(ibeg,jbeg,kbeg, :))
00786             chmax (ibeg,  jv,kv) = maxval (corners(ibeg,jbeg,kbeg, :))
00787          endif
00788 !
00789          j = jend-1; jv = jend
00790          k = kbeg;   kv = kbeg-1 ! ; kc = kbeg
00791 !
00792          if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00793              chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00794             chmin (ibeg,  jv,kv) = array (iend,jbeg,kbeg)
00795             chmax (ibeg-1,jv,kv) = array (iend,jbeg,kbeg)
00796 !
00797             chmin (ibeg-1,jv,kv) = minval (corners(iend,jbeg,kbeg, :))
00798             chmax (ibeg,  jv,kv) = maxval (corners(iend,jbeg,kbeg, :))
00799          endif
00800 !
00801          j = jbeg;   jv = jbeg-1
00802          k = kend-1; kv = kend   ! ; kc = kend
00803 !
00804          if (chmin (i,jend,k) >= minval (corners(i,jend,k, :)) .and. &
00805              chmax (i,jend,k) <= maxval (corners(i,jend,k, :))) then
00806             chmin (iend,  jv,kv) = array (ibeg,jend,kend)
00807             chmax (iend-1,jv,kv) = array (ibeg,jend,kend)
00808 !
00809             chmin (iend-1,jv,kv) = minval (corners(ibeg,jend,kend, :))
00810             chmax (iend,  jv,kv) = maxval (corners(ibeg,jend,kend, :))
00811          endif
00812 !
00813          j = jend-1; jv = jend
00814          k = kend-1; kv = kend ! ; kc = kend
00815 !
00816          if (chmin (i,jend,k) >= minval (corners(i,jend,k, :)) .and. &
00817              chmax (i,jend,k) <= maxval (corners(i,jend,k, :))) then
00818             chmin (iend,  jv,kv) = array (iend,jend,kend)
00819             chmax (iend-1,jv,kv) = array (iend,jend,kend)
00820 !
00821             chmin (iend-1,jv,kv) = minval (corners(iend,jend,kend, :))
00822             chmax (iend,  jv,kv) = maxval (corners(iend,jend,kend, :))
00823          endif
00824 !
00825 !        Lower and upper Face
00826 !
00827             do k = kbeg, kend-1
00828 !cdir vector
00829                do j = jbeg, jend-1
00830                if (chmin (ibeg,j,k) >= minval (corners(ibeg,j,k, :)) .and. &
00831                    chmax (ibeg,j,k) <= maxval (corners(ibeg,j,k, :))) then
00832                   chmin (ibeg,  j, k) = (chmin(ibeg,j,k) + chmax(ibeg,j,k)) * 0.5
00833                   chmax (ibeg-1,j, k) = chmin (ibeg,j,k)
00834 !
00835                   chmin (ibeg-1,j, k) = lchmin (j,k)
00836                   chmax (ibeg,  j, k) = lchmax (j,k)
00837                endif
00838                end do
00839             end do
00840 !
00841          Deallocate (lchmin, lchmax)
00842 !
00843       else
00844 !
00845 !=======================================================================
00846 !     Regular 3d-array
00847 !=======================================================================
00848 !
00849 !     Compute minimum extension of method cell
00850 !
00851          do k = kbeg, kend-1
00852             do j = jbeg, jend-1
00853 !cdir vector
00854                do i = ibeg, iend-1
00855                chmin (i, j, k) = &
00856                  min (array(i,j  ,k  ), array (i+1,j  ,k  ), &
00857                       array(i,j+1,k  ), array (i+1,j+1,k  ), &
00858                       array(i,j  ,k+1), array (i+1,j  ,k+1), &
00859                       array(i,j+1,k+1), array (i+1,j+1,k+1))
00860                enddo
00861             enddo
00862          enddo
00863 !
00864 !     Compute maximum extension of method cell
00865 !
00866          do k = kbeg, kend-1
00867             do j = jbeg, jend-1
00868 !cdir vector
00869                do i = ibeg, iend-1
00870                chmax (i, j, k) = &
00871                  max (array(i,j  ,k  ), array (i+1,j  ,k  ), &
00872                       array(i,j+1,k  ), array (i+1,j+1,k  ), &
00873                       array(i,j  ,k+1), array (i+1,j  ,k+1), &
00874                       array(i,j+1,k+1), array (i+1,j+1,k+1))
00875                enddo
00876             enddo
00877          enddo
00878 !
00879 !     Compute mid point of method cell
00880 !
00881          do k = kbeg, kend-1
00882             do j = jbeg, jend-1
00883 !cdir vector
00884                do i = ibeg, iend-1
00885                midp (i, j, k)  = &
00886                     (array(i,j  ,k  ) + array (i+1,j  ,k  ) + &
00887                      array(i,j+1,k  ) + array (i+1,j+1,k  ) + &
00888                      array(i,j  ,k+1) + array (i+1,j  ,k+1) + &
00889                      array(i,j+1,k+1) + array (i+1,j+1,k+1)) * r_nbr
00890                enddo
00891             enddo
00892          enddo
00893 !
00894 !  ...... Set bounding box of virtual method cells
00895 !
00896 !  ......... for Face 1 
00897 !            virtual cells have k-index kbeg-1
00898 !            corresponding corner index is kbeg
00899 !
00900             do j = jbeg, jend-1
00901 !cdir vector
00902                do i = ibeg, iend-1
00903                chmin (i, j, kbeg-1) = &
00904                  min (minval (corners(i,   j,   kbeg, :)), &
00905                       minval (corners(i+1, j,   kbeg, :)), &
00906                       minval (corners(i,   j+1, kbeg, :)), &
00907                       minval (corners(i+1, j+1, kbeg, :)))
00908                enddo
00909 !
00910 !cdir vector
00911                do i = ibeg, iend-1
00912                chmax (i, j, kbeg-1) = &
00913                  max (maxval (corners(i,   j,   kbeg, :)), &
00914                       maxval (corners(i+1, j,   kbeg, :)), &
00915                       maxval (corners(i,   j+1, kbeg, :)), &
00916                       maxval (corners(i+1, j+1, kbeg, :)))
00917                enddo
00918 !
00919 !cdir vector
00920                do i = ibeg, iend-1
00921                if (chmin (i,j,kbeg-1) <= chmin (i,j,kbeg) .and. &
00922                    chmax (i,j,kbeg-1) >= chmax (i,j,kbeg)) then
00923 !                  chmin (i,j,kbeg-1) =  chmin (i,j,kbeg)
00924 !                  chmax (i,j,kbeg-1) =  chmax (i,j,kbeg)
00925                    chmin (i,j,kbeg-1) = &
00926                      min (array(i,j  ,kbeg), array(i+1,j  ,kbeg), &
00927                           array(i,j+1,kbeg), array(i+1,j+1,kbeg))
00928                    chmax (i,j,kbeg-1) = &
00929                      max (array(i,j  ,kbeg), array(i+1,j  ,kbeg), &
00930                           array(i,j+1,kbeg), array(i+1,j+1,kbeg))
00931                else if (chmin (i,j,kbeg-1) <= chmin (i,j,kbeg)) then
00932                    chmax (i,j,kbeg-1)  = chmin (i,j,kbeg)
00933                else
00934                    chmin (i,j,kbeg-1)  = chmax (i,j,kbeg)
00935                endif
00936                enddo
00937 !
00938             enddo
00939 !
00940 !  ......... for Face 2
00941 !            virtual cells have k-index kend
00942 !            corresponding corner index is kend
00943 !
00944             do j = jbeg, jend-1
00945 !cdir vector
00946                do i = ibeg, iend-1
00947                chmin (i, j, kend) = &
00948                  min (minval (corners(i,   j,   kend, :)), &
00949                       minval (corners(i+1, j,   kend, :)), &
00950                       minval (corners(i,   j+1, kend, :)), &
00951                       minval (corners(i+1, j+1, kend, :)))
00952                enddo
00953 !
00954 !cdir vector
00955                do i = ibeg, iend-1
00956                chmax (i, j, kend) = &
00957                  max (maxval (corners(i,   j,   kend, :)), &
00958                       maxval (corners(i+1, j,   kend, :)), &
00959                       maxval (corners(i,   j+1, kend, :)), &
00960                       maxval (corners(i+1, j+1, kend, :)))
00961                enddo
00962 !
00963 !cdir vector
00964                do i = ibeg, iend-1
00965                if (chmin (i,j,kend) <= chmin (i,j,kend-1) .and. &
00966                    chmax (i,j,kend) >= chmax (i,j,kend-1)) then
00967 !                  chmin (i,j,kend) =  chmin (i,j,kend-1)
00968 !                  chmax (i,j,kend) =  chmax (i,j,kend-1)
00969 
00970                    chmin (i,j,kend) = &
00971                      min (array(i,j  ,kend), array(i+1,j  ,kend), &
00972                           array(i,j+1,kend), array(i+1,j+1,kend))
00973                    chmax (i,j,kend) = &
00974                      max (array(i,j  ,kend), array(i+1,j  ,kend), &
00975                           array(i,j+1,kend), array(i+1,j+1,kend))
00976 
00977                else if (chmin (i,j,kend) <= chmin (i,j,kend-1)) then
00978                    chmax (i,j,kend)  = chmin (i,j,kend-1)
00979                else
00980                    chmin (i,j,kend)  = chmax (i,j,kend-1)
00981                endif
00982                enddo
00983 !
00984             enddo
00985 !
00986 !  ......... for Face 3
00987 !            virtual cells have j-index jbeg-1
00988 !            corresponding corner index is jbeg
00989 !
00990             do k = kbeg, kend-1
00991 !cdir vector
00992                do i = ibeg, iend-1
00993                chmin (i, jbeg-1, k) = &
00994                  min (minval (corners(i,   jbeg, k,   :)), &
00995                       minval (corners(i+1, jbeg, k,   :)), &
00996                       minval (corners(i,   jbeg, k+1, :)), &
00997                       minval (corners(i+1, jbeg, k+1, :)))
00998                enddo
00999 !
01000 !cdir vector
01001                do i = ibeg, iend-1
01002                chmax (i, jbeg-1, k) = &
01003                  max (maxval (corners(i,   jbeg, k,   :)), &
01004                       maxval (corners(i+1, jbeg, k,   :)), &
01005                       maxval (corners(i,   jbeg, k+1, :)), &
01006                       maxval (corners(i+1, jbeg, k+1, :)))
01007                enddo
01008 !
01009 !cdir vector
01010                do i = ibeg, iend-1
01011                if (chmin (i,jbeg-1,k) <= chmin (i,jbeg,k) .and. &
01012                    chmax (i,jbeg-1,k) >= chmax (i,jbeg,k)) then
01013 !                  chmin (i,jbeg-1,k) =  chmin (i,jbeg,k)
01014 !                  chmax (i,jbeg-1,k) =  chmax (i,jbeg,k)
01015 
01016                    chmin (i,jbeg-1,k) = &
01017                      min (array(i,jbeg,k  ), array(i+1,jbeg,k  ), &
01018                           array(i,jbeg,k+1), array(i+1,jbeg,k+1))
01019                    chmax (i,jbeg-1,k) = &
01020                      max (array(i,jbeg,k  ), array(i+1,jbeg,k  ), &
01021                           array(i,jbeg,k+1), array(i+1,jbeg,k+1))
01022 
01023                else if (chmin (i,jbeg-1,k) <= chmin (i,jbeg,k)) then
01024                    chmax (i,jbeg-1,k)  = chmin (i,jbeg,k)
01025                else
01026                    chmin (i,jbeg-1,k)  = chmax (i,jbeg,k)
01027                endif
01028                enddo
01029 !
01030             enddo
01031 !
01032 !  ......... for Face 4
01033 !            virtual cells have j-index jend
01034 !            corresponding corner index is jend
01035 !
01036             do k = kbeg, kend-1
01037 !cdir vector
01038                do i = ibeg, iend-1
01039                chmin (i, jend, k) = &
01040                  min (minval (corners(i,   jend, k,   :)), &
01041                       minval (corners(i+1, jend, k,   :)), &
01042                       minval (corners(i,   jend, k+1, :)), &
01043                       minval (corners(i+1, jend, k+1, :)))
01044                enddo
01045 !
01046 !cdir vector
01047                do i = ibeg, iend-1
01048                chmax (i, jend, k) = &
01049                  max (maxval (corners(i,   jend, k,   :)), &
01050                       maxval (corners(i+1, jend, k,   :)), &
01051                       maxval (corners(i,   jend, k+1, :)), &
01052                       maxval (corners(i+1, jend, k+1, :)))
01053                enddo
01054 !
01055 ! eigentlich will ich die Zelle chmin, chmax nehmen und
01056 ! alles abziehen, was schon vorher in irgendwelchen Zellen ist
01057 ! Fuer global Gitter mit Offset braucht man wirklich die
01058 ! die Methoden Info der benachbarten Bloecke.
01059 !
01060 !cdir vector
01061                do i = ibeg, iend-1
01062                if (chmin (i,jend,k) <= chmin (i,jend-1,k) .and. &
01063                    chmax (i,jend,k) >= chmax (i,jend-1,k)) then
01064 !                  chmin (i,jend,k) =  chmin (i,jend-1,k)
01065 !                  chmax (i,jend,k) =  chmax (i,jend-1,k)
01066 
01067                    chmin (i,jend,k) = &
01068                      min (array(i,jend,k  ), array(i+1,jend,k  ), &
01069                           array(i,jend,k+1), array(i+1,jend,k+1))
01070                    chmax (i,jend,k) = &
01071                      max (array(i,jend,k  ), array(i+1,jend,k  ), &
01072                           array(i,jend,k+1), array(i+1,jend,k+1))
01073 
01074                else if (chmin (i,jend,k) <= chmin (i,jend-1,k)) then
01075                    chmax (i,jend,k)  = chmin (i,jend-1,k)
01076                else
01077                    chmin (i,jend,k)  = chmax (i,jend-1,k)
01078                endif
01079                enddo
01080 
01081             enddo
01082 !
01083 !  ......... for Face 5
01084 !            virtual cells have i-index ibeg-1
01085 !            corresponding corner index is ibeg
01086 !
01087             do k = kbeg, kend-1
01088 !cdir vector
01089                do j = jbeg, jend-1
01090                chmin (ibeg-1, j, k) = &
01091                  min (minval (corners(ibeg, j,     k,   :)), &
01092                       minval (corners(ibeg, j+1,   k,   :)), &
01093                       minval (corners(ibeg, j,     k+1, :)), &
01094                       minval (corners(ibeg, j+1,   k+1, :)))
01095                enddo
01096 !
01097 !cdir vector
01098                do j = jbeg, jend-1
01099                chmax (ibeg-1, j, k) = &
01100                  max (maxval (corners(ibeg, j,     k,   :)), &
01101                       maxval (corners(ibeg, j+1,   k,   :)), &
01102                       maxval (corners(ibeg, j,     k+1, :)), &
01103                       maxval (corners(ibeg, j+1,   k+1, :)))
01104                enddo
01105 !
01106 !cdir vector
01107                do j = jbeg, jend-1
01108                if (chmin (ibeg-1,j,k) <= chmin (ibeg,j,k) .and. &
01109                    chmax (ibeg-1,j,k) >= chmax (ibeg,j,k)) then
01110 !                  chmin (ibeg-1,j,k) =  chmin (ibeg,j,k)
01111 !                  chmax (ibeg-1,j,k) =  chmax (ibeg,j,k)
01112 
01113                    chmin (ibeg-1,j,k) = &
01114                      min (array(ibeg,j,k  ), array(ibeg,j+1,k  ), &
01115                           array(ibeg,j,k+1), array(ibeg,j+1,k+1))
01116                    chmax (ibeg-1,j,k) = &
01117                      max (array(ibeg,j,k  ), array(ibeg,j+1,k  ), &
01118                           array(ibeg,j,k+1), array(ibeg,j+1,k+1))
01119 
01120                else if (chmin (ibeg-1,j,k) <= chmin (ibeg,j,k)) then
01121                    chmax (ibeg-1,j,k)  = chmin (ibeg,j,k)
01122                else
01123                    chmin (ibeg-1,j,k)  = chmax (ibeg,j,k)
01124                endif
01125                enddo
01126 !
01127             enddo
01128 !
01129 !  ......... for Face 6
01130 !            virtual cells have i-index iend
01131 !            corresponding corner index is iend
01132 !
01133             do k = kbeg, kend-1
01134 !cdir vector
01135                do j = jbeg, jend-1
01136                chmin (iend, j, k) = &
01137                  min (minval (corners(iend, j,     k,   :)), &
01138                       minval (corners(iend, j+1,   k,   :)), &
01139                       minval (corners(iend, j,     k+1, :)), &
01140                       minval (corners(iend, j+1,   k+1, :)))
01141                enddo
01142 !
01143 !cdir vector
01144                do j = jbeg, jend-1
01145                chmax (iend, j, k) = &
01146                  max (maxval (corners(iend, j,     k,   :)), &
01147                       maxval (corners(iend, j+1,   k,   :)), &
01148                       maxval (corners(iend, j,     k+1, :)), &
01149                       maxval (corners(iend, j+1,   k+1, :)))
01150                enddo
01151 !
01152 !cdir vector
01153                do j = jbeg, jend-1
01154                if (chmin (iend,j,k) <= chmin (iend-1,j,k) .and. &
01155                    chmax (iend,j,k) >= chmax (iend-1,j,k)) then
01156 !                  chmin (iend,j,k) =  chmin (iend-1,j,k)
01157 !                  chmax (iend,j,k) =  chmax (iend-1,j,k)
01158 
01159                    chmin (iend,j,k) = &
01160                      min (array(iend,j,k  ), array(iend,j+1,k  ), &
01161                           array(iend,j,k+1), array(iend,j+1,k+1))
01162                    chmax (iend,j,k) = &
01163                      max (array(iend,j,k  ), array(iend,j+1,k  ), &
01164                           array(iend,j,k+1), array(iend,j+1,k+1))
01165 
01166                else if (chmin (iend,j,k) <= chmin (iend-1,j,k)) then
01167                    chmax (iend,j,k)  = chmin (iend-1,j,k)
01168                else
01169                    chmin (iend,j,k)  = chmax (iend-1,j,k)
01170                endif
01171                enddo
01172             enddo
01173 !
01174 !  ......... for Edges of Faces
01175 !
01176 !        Edge of faces 1 and 3 (kbeg   and jbeg)
01177 !        Edge of faces 1 and 4 (kbeg   and jend-1)
01178 !        Edge of faces 2 and 3 (kend-1 and jbeg)
01179 !        Edge of faces 2 and 4 (kend-1 and jend-1)
01180 !
01181 ! Um weniger Code zu haben: ??? Performance ???
01182 ! kvt = Index of       virtual method cell - in k-direction
01183 ! kc  = Index of corresponding corner cell - in k-direction
01184 ! kin = Index of last interior method cell - in k-direction
01185 ! jvt = Index of       virtual method cell - in j-direction
01186 ! jc  = Index of corresponding corner cell - in j-direction
01187 ! jin = Index of last interior method cell - in j-direction
01188 !
01189          kin (1) = kbeg;   kvt (1) = kbeg - 1; kc (1) = kbeg
01190          jin (1) = jbeg;   jvt (1) = jbeg - 1; jc (1) = jbeg
01191 !
01192          kin (2) = kbeg;   kvt (2) = kbeg - 1; kc (2) = kbeg
01193          jin (2) = jend-1; jvt (2) = jend;     jc (2) = jend
01194 !
01195          kin (3) = kend-1; kvt (3) = kend    ; kc (3) = kend
01196          jin (3) = jbeg;   jvt (3) = jbeg - 1; jc (3) = jbeg
01197 !
01198          kin (4) = kend-1; kvt (4) = kend    ; kc (4) = kend
01199          jin (4) = jend-1; jvt (4) = jend;     jc (4) = jend
01200 !
01201             do n = 1, 4
01202 !cdir vector
01203                do i = ibeg, iend-1
01204                chmin (i, jvt(n), kvt(n)) = &
01205                  min (minval (corners(i,   jc(n), kc(n), :)), &
01206                       minval (corners(i+1, jc(n), kc(n), :)))
01207                enddo
01208 !
01209 !cdir vector
01210                do i = ibeg, iend-1
01211                chmax (i, jvt(n), kvt(n)) = &
01212                  max (maxval (corners(i,   jc(n), kc(n), :)), &
01213                       maxval (corners(i+1, jc(n), kc(n), :)))
01214                enddo
01215 !
01216 !cdir vector
01217                do i = ibeg, iend-1
01218                if (chmin (i, jvt(n),kvt(n)) <= chmin (i,jin(n),kin(n)) .and. &
01219                    chmax (i, jvt(n),kvt(n)) >= chmax (i,jin(n),kin(n))) then
01220 
01221                    chmin (i, jvt(n),kvt(n))  = chmin (i,jin(n),kin(n))
01222                    chmax (i, jvt(n),kvt(n))  = chmax (i,jin(n),kin(n))
01223 
01224                else if (chmin (i, jvt(n),kvt(n)) <= chmin (i,jin(n),kin(n))) then
01225                    chmax (i, jvt(n),kvt(n))  = chmin (i,jin(n),kin(n))
01226                else
01227                    chmin (i, jvt(n),kvt(n))  = chmax (i,jin(n),kin(n))
01228                endif
01229                enddo
01230             end do
01231 !
01232 !        Edge of faces 1 and 5 (kbeg   and ibeg)
01233 !        Edge of faces 1 and 6 (kbeg   and iend-1)
01234 !        Edge of faces 2 and 5 (kend-1 and ibeg)
01235 !        Edge of faces 2 and 6 (kend-1 and iend-1)
01236 !
01237 !        kin (1) = kbeg;   kvt (1) = kbeg - 1; kc (1) = kbeg
01238          iin (1) = ibeg;   ivt (1) = ibeg - 1; ic (1) = ibeg
01239 !
01240 !        kin (2) = kbeg;   kvt (2) = kbeg - 1; kc (2) = kbeg
01241          iin (2) = iend-1; ivt (2) = iend    ; ic (2) = iend
01242 !
01243 !        kin (3) = kend-1; kvt (3) = kend    ; kc (3) = kend
01244          iin (3) = ibeg;   ivt (3) = ibeg - 1; ic (3) = ibeg
01245 !
01246 !        kin (4) = kend-1; kvt (4) = kend    ; kc (4) = kend
01247          iin (4) = iend-1; ivt (4) = iend    ; ic (4) = iend
01248 !
01249             do n = 1, 4
01250 !
01251 !cdir vector
01252                do j = jbeg, jend-1
01253                chmin (ivt(n), j, kvt(n)) = &
01254                  min (minval (corners(ic(n), j,   kc(n), :)), &
01255                       minval (corners(ic(n), j+1, kc(n), :)))
01256                enddo
01257 !
01258 !cdir vector
01259                do j = jbeg, jend-1
01260                chmax (ivt(n), j, kvt(n)) = &
01261                  max (maxval (corners(ic(n), j,   kc(n), :)), &
01262                       maxval (corners(ic(n), j+1, kc(n), :)))
01263                enddo
01264 !
01265 !cdir vector
01266                do j = jbeg, jend-1
01267                if (chmin (ivt(n), j, kvt(n)) <= chmin (iin(n),j,kin(n)) .and. &
01268                    chmax (ivt(n), j, kvt(n)) >= chmax (iin(n),j,kin(n))) then
01269 
01270                    chmin (ivt(n), j, kvt(n))  = chmin (iin(n),j,kin(n))
01271                    chmax (ivt(n), j, kvt(n))  = chmax (iin(n),j,kin(n))
01272 
01273                else if (chmin (ivt(n), j, kvt(n)) <= chmin (iin(n),j,kin(n))) then
01274                    chmax (ivt(n), j, kvt(n))  = chmin (iin(n),j,kin(n))
01275                else
01276                    chmin (ivt(n), j, kvt(n))  = chmax (iin(n),j,kin(n))
01277                endif
01278                enddo
01279             end do
01280 !
01281 !        Edge of faces 2 and 5 (jbeg   and ibeg)
01282 !        Edge of faces 2 and 6 (jbeg   and iend-1)
01283 !        Edge of faces 3 and 5 (jend-1 and ibeg)
01284 !        Edge of faces 3 and 6 (jend-1 and iend-1)
01285 !
01286          jin (1) = jbeg;   jvt (1) = jbeg - 1; jc (1) = jbeg
01287 !        iin (1) = ibeg;   ivt (1) = ibeg - 1; ic (1) = ibeg
01288 !
01289          jin (2) = jbeg;   jvt (2) = jbeg - 1; jc (2) = jbeg
01290 !        iin (2) = iend-1; ivt (2) = iend    ; ic (2) = iend
01291 !
01292          jin (3) = jend-1; jvt (3) = jend;     jc (3) = jend
01293 !        iin (3) = ibeg;   ivt (3) = ibeg - 1; ic (3) = ibeg
01294 !
01295          jin (4) = jend-1; jvt (4) = jend;     jc (4) = jend
01296 !        iin (4) = iend-1; ivt (4) = iend    ; ic (4) = iend
01297 !
01298             do n = 1, 4
01299 !
01300 !cdir vector
01301                do k = kbeg, kend-1
01302                chmin (ivt(n), jvt(n), k) = &
01303                  min (minval (corners(ic(n), jc(n), k  , :)), &
01304                       minval (corners(ic(n), jc(n), k+1, :)))
01305                enddo
01306 !
01307 !cdir vector
01308                do k = kbeg, kend-1
01309                chmax (ivt(n), jvt(n), k) = &
01310                  max (maxval (corners(ic(n), jc(n), k  , :)), &
01311                       maxval (corners(ic(n), jc(n), k+1, :)))
01312                enddo
01313 !
01314 !cdir vector
01315                do k = kbeg, kend-1
01316                if (chmin (ivt(n),jvt(n), k) <= chmin (iin(n),jin(n),k) .and. &
01317                    chmax (ivt(n),jvt(n), k) >= chmax (iin(n),jin(n),k)) then
01318 
01319                    chmin (ivt(n),jvt(n), k)  = chmin (iin(n),jin(n),k)
01320                    chmax (ivt(n),jvt(n), k)  = chmax (iin(n),jin(n),k)
01321 
01322                else if (chmin (ivt(n),jvt(n), k) <= chmin (iin(n),jin(n),k)) then
01323                    chmax (ivt(n),jvt(n), k)  = chmin (iin(n),jin(n),k)
01324                else
01325                    chmin (ivt(n),jvt(n), k)  = chmax (iin(n),jin(n),k)
01326                endif
01327                enddo
01328             end do
01329 !
01330 !  ......... for Corners of edges
01331 !
01332 ! kvt = Index of virtual method cell - in k-direction
01333 ! kin = Last interior index          - in k-direction
01334 ! jvt = Index of virtual method cell - in j-direction
01335 ! jin = Last interior index          - in j-direction
01336 !
01337          kin (1) = kbeg;   kvt (1) = kbeg - 1; kc (1) = kbeg
01338          jin (1) = jbeg;   jvt (1) = jbeg - 1; jc (1) = jbeg
01339          iin (1) = ibeg;   ivt (1) = ibeg - 1; ic (1) = ibeg
01340 !
01341          kin (2) = kbeg;   kvt (2) = kbeg - 1; kc (2) = kbeg
01342          jin (2) = jbeg;   jvt (2) = jbeg - 1; jc (2) = jbeg
01343          iin (2) = iend-1; ivt (2) = iend    ; ic (2) = iend
01344 !
01345          kin (3) = kbeg;   kvt (3) = kbeg - 1; kc (3) = kbeg
01346          jin (3) = jend-1; jvt (3) = jend;     jc (3) = jend
01347          iin (3) = ibeg;   ivt (3) = ibeg - 1; ic (3) = ibeg
01348 !
01349          kin (4) = kbeg;   kvt (4) = kbeg - 1; kc (4) = kbeg
01350          jin (4) = jend-1; jvt (4) = jend;     jc (4) = jend
01351          iin (4) = iend-1; ivt (4) = iend    ; ic (4) = iend
01352 !
01353          kin (5) = kend-1; kvt (5) = kend    ; kc (5) = kend
01354          jin (5) = jbeg;   jvt (5) = jbeg - 1; jc (5) = jbeg
01355          iin (5) = ibeg;   ivt (5) = ibeg - 1; ic (5) = ibeg
01356 !
01357          kin (6) = kend-1; kvt (6) = kend    ; kc (6) = kend
01358          jin (6) = jbeg;   jvt (6) = jbeg - 1; jc (6) = jbeg
01359          iin (6) = iend-1; ivt (6) = iend    ; ic (6) = iend
01360 !
01361          kin (7) = kend-1; kvt (7) = kend    ; kc (7) = kend
01362          jin (7) = jend-1; jvt (7) = jend;     jc (7) = jend
01363          iin (7) = ibeg;   ivt (7) = ibeg - 1; ic (7) = ibeg
01364 !
01365          kin (n_corners_3d) = kend-1
01366          kvt (n_corners_3d) = kend
01367          kc  (n_corners_3d) = kend
01368          jin (n_corners_3d) = jend-1
01369          jvt (n_corners_3d) = jend
01370          jc  (n_corners_3d) = jend
01371          iin (n_corners_3d) = iend-1
01372          ivt (n_corners_3d) = iend
01373          ic  (n_corners_3d) = iend
01374 !
01375             do n = 1, n_corners_3d
01376             chmin (ivt(n), jvt(n), kvt(n)) = &
01377                minval (corners(ic(n), jc(n), kc(n), :))
01378             chmax (ivt(n), jvt(n), kvt(n)) = &
01379                maxval (corners(ic(n), jc(n), kc(n), :))
01380 !
01381             if (chmin (ivt(n),jvt(n),kvt(n)) <= chmin (iin(n),jin(n),kin(n)) .and. &
01382                 chmax (ivt(n),jvt(n),kvt(n)) >= chmax (iin(n),jin(n),kin(n))) then
01383 
01384                 chmin (ivt(n),jvt(n),kvt(n))  = chmin (iin(n),jin(n),kin(n))
01385                 chmax (ivt(n),jvt(n),kvt(n))  = chmax (iin(n),jin(n),kin(n))
01386 
01387             else if (chmin (ivt(n),jvt(n),kvt(n)) <= chmin (iin(n),jin(n),kin(n))) &
01388             then
01389                 chmax (ivt(n),jvt(n),kvt(n))  = chmin (iin(n),jin(n),kin(n))
01390             else
01391                 chmin (ivt(n),jvt(n),kvt(n))  = chmax (iin(n),jin(n),kin(n))
01392             endif
01393             end do
01394 !
01395 !  ... Set mid point of virtual method cells
01396 !
01397          k = kbeg - 1
01398             do j = jbeg-1, jend
01399 !cdir vector
01400                do i = ibeg-1, iend
01401                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01402                enddo
01403             enddo
01404 !
01405          k = kend
01406             do j = jbeg-1, jend
01407 !cdir vector
01408                do i = ibeg-1, iend
01409                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01410                enddo
01411             enddo
01412 !
01413 !
01414 !
01415          j = jbeg - 1
01416             do k = kbeg, kend-1
01417 !cdir vector
01418                do i = ibeg-1, iend
01419                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01420                enddo
01421             enddo
01422 !
01423          j = jend
01424             do k = kbeg, kend-1
01425 !cdir vector
01426                do i = ibeg-1, iend
01427                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01428                enddo
01429             enddo
01430 !
01431 !
01432 !
01433          i = ibeg - 1
01434             do k = kbeg, kend-1
01435 !cdir vector
01436                do j = jbeg, jend-1
01437                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01438                enddo
01439             enddo
01440 !
01441          i = iend
01442             do k = kbeg, kend-1
01443 !cdir vector
01444                do j = jbeg, jend-1
01445                midp  (i, j, k) = (chmin(i,j,k) + chmax (i,j,k)) * 0.5
01446                enddo
01447             enddo
01448 
01449       endif
01450 !
01451 !===> All done
01452 !
01453 #ifdef VERBOSE
01454       print 9980, trim(ch_id), ierror
01455 
01456       call psmile_flushstd
01457 #endif /* VERBOSE */
01458 !
01459 !  Formats:
01460 !
01461 9990 format (1x, a, ': psmile_bbcells_3d_dble')
01462 9980 format (1x, a, ': psmile_bbcells_3d_dble: eof ierror =', i3)
01463 !
01464       end subroutine PSMILe_Bbcells_3d_dble
Generated on 18 Mar 2011 for Oasis4 by  doxygen 1.6.1