psmile_get_cyclic_dir_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_Get_cyclic_dir_3d_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_get_cyclic_dir_3d_dble (               &
00012                            chmin, chmin2, chmin3,              &
00013                            chmax, chmax2, chmax3, levdim,      &
00014                            corners, corners2, corners3,        &
00015                            corner_shape, nbr_corners,          &
00016                            grid_valid_shape,                   &
00017                            len_cyclic, rtol, index, cyclic, ierror)
00018 !
00019 ! !USES:
00020 !
00021       use PRISM_constants
00022 !
00023       use PSMILe, dummy_interface => PSMILe_Get_cyclic_dir_3d_dble
00024 
00025       implicit none
00026 !
00027 ! !INPUT PARAMETERS:
00028 !
00029 
00030       Integer, Intent (In)            :: levdim (ndim_3d)
00031 
00032 !     Dimension of chmin, chmax
00033 
00034       Double Precision, Intent (In)   :: chmin  (0:levdim(1), 
00035                                                  0:levdim(2), 
00036                                                  0:levdim(3))
00037       Double Precision, Intent (In)   :: chmin2 (0:levdim(1), 
00038                                                  0:levdim(2), 
00039                                                  0:levdim(3))
00040       Double Precision, Intent (In)   :: chmin3 (0:levdim(1), 
00041                                                  0:levdim(2), 
00042                                                  0:levdim(3))
00043 
00044 !     Minimum of grid coordinates per grid cell
00045 
00046       Double Precision, Intent (In)   :: chmax  (0:levdim(1), 
00047                                                  0:levdim(2), 
00048                                                  0:levdim(3))
00049       Double Precision, Intent (In)   :: chmax2 (0:levdim(1), 
00050                                                  0:levdim(2), 
00051                                                  0:levdim(3))
00052       Double Precision, Intent (In)   :: chmax3 (0:levdim(1), 
00053                                                  0:levdim(2), 
00054                                                  0:levdim(3))
00055 
00056 !     Maximum of grid coordinates per grid cell
00057 
00058 !
00059       Integer, Intent (In)            :: corner_shape (2, ndim_3d)
00060 
00061 !     Dimension of corner coordinates
00062 !
00063       Integer, Intent (In)            :: nbr_corners
00064 !
00065 !     Number of corners per corner cell
00066 
00067       Double Precision, Intent (In)   :: corners(corner_shape(1,1): 
00068                                                   corner_shape(2,1), 
00069                                                   corner_shape(1,2): 
00070                                                   corner_shape(2,2), 
00071                                                   corner_shape(1,3): 
00072                                                   corner_shape(2,3), 
00073                                                   nbr_corners)
00074       Double Precision, Intent (In)   :: corners2(corner_shape(1,1): 
00075                                                   corner_shape(2,1), 
00076                                                   corner_shape(1,2): 
00077                                                   corner_shape(2,2), 
00078                                                   corner_shape(1,3): 
00079                                                   corner_shape(2,3), 
00080                                                   nbr_corners)
00081       Double Precision, Intent (In)   :: corners3(corner_shape(1,1): 
00082                                                   corner_shape(2,1), 
00083                                                   corner_shape(1,2): 
00084                                                   corner_shape(2,2), 
00085                                                   corner_shape(1,3): 
00086                                                   corner_shape(2,3), 
00087                                                   nbr_corners)
00088 
00089 !     Corner-Coordinates to be controlled
00090 !
00091 
00092       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00093 
00094 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00095 
00096       Double Precision, Intent (In)   :: len_cyclic
00097 
00098 !     Cycle-length of coordinate 
00099 
00100       Double Precision, Intent (In)   :: rtol
00101 
00102 !     Tolerance
00103 
00104       Integer, Intent (In)            :: index
00105 
00106 !     Specifies the direction (1(i), 2(j), 3(k)) to be controlled
00107 !
00108 ! !OUTPUT PARAMETERS:
00109 !
00110       Logical, Intent (Out)           :: cyclic
00111 !
00112 !     Is this corner coordinate a cyclic coordinate in the direction ?
00113 !     cyclic = .true.   : Coordinate is     a cyclic coordinate
00114 !     cyclic = .false.. : Coordinate is not a cyclic coordinate
00115 !
00116       Integer, Intent (Out)           :: ierror
00117 
00118 !     Returns the error code of PSMILE_Get_cyclic_dir_3d_dble;
00119 !             ierror = 0 : No error
00120 !             ierror > 0 : Severe error
00121 !
00122 ! !LOCAL VARIABLES
00123 !
00124 !     ... for comparison
00125 !
00126       Integer                         :: i, j, k
00127       Double Precision                :: dmin
00128       Double Precision, Allocatable   :: face1_min (:, :), face1_max (:, :)
00129       Double Precision, Allocatable   :: face2_min (:, :), face2_max (:, :)
00130       Double Precision, Allocatable   :: min1 (:, :), max1 (:, :)
00131 !
00132 !     ... for error handling
00133 !
00134       Integer, parameter              :: nerrp = 1
00135       Integer                         :: ierrp (nerrp)
00136 !
00137 !
00138 ! !DESCRIPTION:
00139 !
00140 ! Subroutine "PSMILe_Get_cyclic_dir_3d_dble" determines whether
00141 ! the current block is cyclic in the given logical direction "index".
00142 !
00143 ! It's assumed that the coordinate to be controlled for the cylic index
00144 ! is stored in "chmin", "chmax" and "corners".
00145 ! The 2 other coordinates "chmin2", "chmin3", ... are used to control
00146 ! whether the cells are fitting.
00147 !
00148 !
00149 ! !REVISION HISTORY:
00150 !
00151 !   Date      Programmer   Description
00152 ! ----------  ----------   -----------
00153 ! 03.07.21    H. Ritzdorf  created
00154 !
00155 !EOP
00156 !----------------------------------------------------------------------
00157 !
00158 !  $Id: psmile_get_cyclic_dir_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $
00159 !  $Author: valcke $
00160 !
00161    Character(len=len_cvs_string), save :: mycvs = 
00162        '$Id: psmile_get_cyclic_dir_3d_dble.F90 2325 2010-04-21 15:00:07Z valcke $'
00163 !
00164 !----------------------------------------------------------------------
00165 !
00166 !  Initialization
00167 !
00168 #ifdef VERBOSE
00169       print 9990, trim(ch_id), index
00170 
00171       call psmile_flushstd
00172 #endif /* VERBOSE */
00173 !
00174 #ifdef PRISM_ASSERTION
00175 #endif
00176 !
00177       ierror = 0
00178       cyclic = .false.
00179 !
00180 !-----------------------------------------------------------------------
00181 !     This is a weak control. But without definition of a
00182 !     volume type and a order of points it's not really possible 
00183 !     to determine wheher index i is really cyclic.
00184 !-----------------------------------------------------------------------
00185 !-----------------------------------------------------------------------------
00186 !     Control logical index 1
00187 !-----------------------------------------------------------------------------
00188 !
00189       if (index == 1) then
00190          dmin = minval ( chmax (levdim(1), :, :) - chmin (0, :, :) )
00191          if (dmin < len_cyclic) return
00192 !
00193             do k = 0, levdim (3)
00194 !cdir vector
00195                do j = 0, levdim (2)
00196                if (chmax2(levdim(1), j, k) < chmin2 (0,         j, k) .or. &
00197                    chmax2(0,         j, k) < chmin2 (levdim(1), j, k)) return
00198                end do
00199             end do
00200 !
00201             do k = 0, levdim (3)
00202 !cdir vector
00203                do j = 0, levdim (2)
00204                if (chmax3(levdim(1), j, k) < chmin3 (0,         j, k) .or. &
00205                    chmax3(0,         j, k) < chmin3 (levdim(1), j, k)) return
00206                end do
00207             end do
00208 !
00209 !===> ... transform possible cyclic index
00210 !
00211             allocate (face1_min (0:levdim(2), 0:levdim(3)), stat = ierror)
00212             if ( ierror /= 0 ) then
00213                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00214                ierror = PRISM_Error_Alloc
00215                call psmile_error ( ierror, 'face1_min', &
00216                                    ierrp, 1, __FILE__, __LINE__ )
00217                return
00218             endif
00219 
00220             allocate (face1_max (0:levdim(2), 0:levdim(3)), stat = ierror)
00221             if ( ierror /= 0 ) then
00222                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00223                ierror = PRISM_Error_Alloc
00224                call psmile_error ( ierror, 'face1_max', &
00225                                    ierrp, 1, __FILE__, __LINE__ )
00226                return
00227             endif
00228 
00229             allocate (face2_min (0:levdim(2), 0:levdim(3)), stat = ierror)
00230             if ( ierror /= 0 ) then
00231                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00232                ierror = PRISM_Error_Alloc
00233                call psmile_error ( ierror, 'face2_min', &
00234                                    ierrp, 1, __FILE__, __LINE__ )
00235                return
00236             endif
00237 
00238             allocate (face2_max (0:levdim(2), 0:levdim(3)), stat = ierror)
00239             if ( ierror /= 0 ) then
00240                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00241                ierror = PRISM_Error_Alloc
00242                call psmile_error ( ierror, 'face2_max', &
00243                                    ierrp, 1, __FILE__, __LINE__ )
00244                return
00245             endif
00246 
00247             allocate (min1 (0:levdim(2), 0:levdim(3)), stat = ierror)
00248             if ( ierror /= 0 ) then
00249                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00250                ierror = PRISM_Error_Alloc
00251                call psmile_error ( ierror, 'min1', &
00252                                    ierrp, 1, __FILE__, __LINE__ )
00253                return
00254             endif
00255 
00256             allocate (max1 (0:levdim(2), 0:levdim(3)), stat = ierror)
00257             if ( ierror /= 0 ) then
00258                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00259                ierror = PRISM_Error_Alloc
00260                call psmile_error ( ierror, 'max1', &
00261                                    ierrp, 1, __FILE__, __LINE__ )
00262                return
00263             endif
00264 !
00265             min1 (:, :) = mod (chmin (0, :, :), len_cyclic)
00266             max1 (:, :) = mod (chmax (0, :, :), len_cyclic)
00267 
00268             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00269             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00270 
00271             face1_min = min (min1, max1) - rtol
00272             face1_max = max (min1, max1)
00273 !
00274             min1 (:, :) = mod (chmin (levdim(1), :, :), len_cyclic)
00275             max1 (:, :) = mod (chmax (levdim(1), :, :), len_cyclic)
00276 
00277             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00278             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00279 
00280             face2_min = min (min1, max1) - rtol
00281             face2_max = max (min1, max1)
00282 !
00283         k1: do k = 0, levdim (3)
00284 !cdir vector
00285                do j = 0, levdim (2)
00286                if (face1_max(j, k) < face2_min (j, k) .or. &
00287                    face2_max(j, k) < face1_min (j, k)) exit k1
00288                end do
00289 
00290             end do k1
00291 !
00292          cyclic = k > levdim (3)
00293 
00294 !-----------------------------------------------------------------------------
00295 !        Control logical index 2
00296 !-----------------------------------------------------------------------------
00297 !
00298       else if (index == 2) then
00299 !
00300          dmin = minval ( chmax (:, levdim(2), :) - chmin (:, 0, :) )
00301          if (dmin < len_cyclic) return
00302 !
00303             do k = 0, levdim (3)
00304 !cdir vector
00305                do i = 0, levdim (1)
00306                if (chmax2(i, levdim(2), k) < chmin2 (i, 0,         k) .or. &
00307                    chmax2(i, 0,         k) < chmin2 (i, levdim(2), k)) return
00308                end do
00309             end do
00310 !
00311             do k = 0, levdim (3)
00312 !cdir vector
00313                do i = 0, levdim (1)
00314                if (chmax3(i, levdim(2), k) < chmin3 (i, 0,         k) .or. &
00315                    chmax3(i, 0,         k) < chmin3 (i, levdim(2), k)) return
00316                end do
00317             end do
00318 !
00319 !===> ... transform possible cyclic index
00320 !
00321             allocate (face1_min (0:levdim(1), 0:levdim(3)), stat = ierror)
00322             if ( ierror /= 0 ) then
00323                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00324                ierror = PRISM_Error_Alloc
00325                call psmile_error ( ierror, 'face1_min', &
00326                                    ierrp, 1, __FILE__, __LINE__ )
00327                return
00328             endif
00329 
00330             allocate (face1_max (0:levdim(1), 0:levdim(3)), stat = ierror)
00331             if ( ierror /= 0 ) then
00332                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00333                ierror = PRISM_Error_Alloc
00334                call psmile_error ( ierror, 'face1_max', &
00335                                    ierrp, 1, __FILE__, __LINE__ )
00336                return
00337             endif
00338 
00339             allocate (face2_min (0:levdim(1), 0:levdim(3)), stat = ierror)
00340             if ( ierror /= 0 ) then
00341                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00342                ierror = PRISM_Error_Alloc
00343                call psmile_error ( ierror, 'face2_min', &
00344                                    ierrp, 1, __FILE__, __LINE__ )
00345                return
00346             endif
00347 
00348             allocate (face2_max (0:levdim(1), 0:levdim(3)), stat = ierror)
00349             if ( ierror /= 0 ) then
00350                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00351                ierror = PRISM_Error_Alloc
00352                call psmile_error ( ierror, 'face2_max', &
00353                                    ierrp, 1, __FILE__, __LINE__ )
00354                return
00355             endif
00356 
00357             allocate (min1 (0:levdim(1), 0:levdim(3)), stat = ierror)
00358             if ( ierror /= 0 ) then
00359                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00360                ierror = PRISM_Error_Alloc
00361                call psmile_error ( ierror, 'min1', &
00362                                    ierrp, 1, __FILE__, __LINE__ )
00363                return
00364             endif
00365 !
00366             allocate (max1 (0:levdim(1), 0:levdim(3)), stat = ierror)
00367             if ( ierror /= 0 ) then
00368                ierrp (1) = (levdim(1)+1) * (levdim(3)+1)
00369                ierror = PRISM_Error_Alloc
00370                call psmile_error ( ierror, 'max1', &
00371                                    ierrp, 1, __FILE__, __LINE__ )
00372                return
00373             endif
00374 !
00375             min1 (:, :) = mod (chmin (:, 0, :), len_cyclic)
00376             max1 (:, :) = mod (chmax (:, 0, :), len_cyclic)
00377 
00378             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00379             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00380 
00381             face1_min = min (min1, max1) - rtol
00382             face1_max = max (min1, max1)
00383 !
00384             min1 (:, :) = mod (chmin (:, levdim(2), :), len_cyclic)
00385             max1 (:, :) = mod (chmax (:, levdim(2), :), len_cyclic)
00386 
00387             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00388             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00389 
00390             face2_min = min (min1, max1) - rtol
00391             face2_max = max (min1, max1)
00392 !
00393         k2: do k = 0, levdim (3)
00394 !cdir vector
00395                do i = 0, levdim (1)
00396                if (face1_max(i, k) < face2_min (i, k) .or. &
00397                    face2_max(i, k) < face1_min (i, k)) exit k2
00398                end do
00399             end do k2
00400 !
00401          cyclic = k > levdim (3)
00402 !
00403       else if (index == 3) then
00404 
00405 !-----------------------------------------------------------------------------
00406 !        Control logical index 3
00407 !-----------------------------------------------------------------------------
00408 !
00409          dmin = minval ( chmax (:, :, levdim(3)) - chmin (:, :, 0) )
00410          if (dmin < len_cyclic) return
00411 !
00412 !===> ... transform possible cyclic index
00413 !
00414             allocate (face1_min (0:levdim(1), 0:levdim(2)), stat = ierror)
00415             if ( ierror /= 0 ) then
00416                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00417                ierror = PRISM_Error_Alloc
00418                call psmile_error ( ierror, 'face1_min', &
00419                                    ierrp, 1, __FILE__, __LINE__ )
00420                return
00421             endif
00422 
00423             allocate (face1_max (0:levdim(1), 0:levdim(2)), stat = ierror)
00424             if ( ierror /= 0 ) then
00425                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00426                ierror = PRISM_Error_Alloc
00427                call psmile_error ( ierror, 'face1_max', &
00428                                    ierrp, 1, __FILE__, __LINE__ )
00429                return
00430             endif
00431 
00432             allocate (face2_min (0:levdim(1), 0:levdim(2)), stat = ierror)
00433             if ( ierror /= 0 ) then
00434                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00435                ierror = PRISM_Error_Alloc
00436                call psmile_error ( ierror, 'face2_min', &
00437                                    ierrp, 1, __FILE__, __LINE__ )
00438                return
00439             endif
00440 
00441             allocate (face2_max (0:levdim(1), 0:levdim(2)), stat = ierror)
00442             if ( ierror /= 0 ) then
00443                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00444                ierror = PRISM_Error_Alloc
00445                call psmile_error ( ierror, 'face2_max', &
00446                                    ierrp, 1, __FILE__, __LINE__ )
00447                return
00448             endif
00449 
00450             allocate (min1 (0:levdim(1), 0:levdim(2)), stat = ierror)
00451             if ( ierror /= 0 ) then
00452                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00453                ierror = PRISM_Error_Alloc
00454                call psmile_error ( ierror, 'min1', &
00455                                    ierrp, 1, __FILE__, __LINE__ )
00456                return
00457             endif
00458 !
00459             allocate (max1 (0:levdim(1), 0:levdim(2)), stat = ierror)
00460             if ( ierror /= 0 ) then
00461                ierrp (1) = (levdim(1)+1) * (levdim(2)+1)
00462                ierror = PRISM_Error_Alloc
00463                call psmile_error ( ierror, 'max1', &
00464                                    ierrp, 1, __FILE__, __LINE__ )
00465                return
00466             endif
00467 !
00468             min1 (:, :) = mod (chmin (:, :, 0), len_cyclic)
00469             max1 (:, :) = mod (chmax (:, :, 0), len_cyclic)
00470 
00471             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00472             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00473 
00474             face1_min = min (min1, max1) - rtol
00475             face1_max = max (min1, max1)
00476 !
00477             min1 (:, :) = mod (chmin (:, :, levdim(3)), len_cyclic)
00478             max1 (:, :) = mod (chmax (:, :, levdim(3)), len_cyclic)
00479 
00480             where (min1 .lt. 0.0) min1 = min1 + len_cyclic
00481             where (max1 .lt. 0.0) max1 = max1 + len_cyclic
00482 
00483             face2_min = min (min1, max1) - rtol
00484             face2_max = max (min1, max1)
00485 !
00486         j3: do j = 0, levdim (2)
00487 !cdir vector
00488                do i = 0, levdim (1)
00489                if (face1_max(i, j) < face2_min (i, j) .or. &
00490                    face2_max(i, j) < face1_min (i, j)) exit j3
00491                end do
00492 
00493             end do j3
00494 !
00495          cyclic = j > levdim (2)
00496 
00497       else
00498 
00499 !-----------------------------------------------------------------------------
00500 !        Internal error
00501 !-----------------------------------------------------------------------------
00502 
00503          ierror = PRISM_Error_Internal
00504          ierrp (1) = index
00505 
00506          call psmile_error ( ierror, 'unsupported index', &
00507                              ierrp, 1, __FILE__, __LINE__ )
00508 
00509       endif
00510 !
00511 !===> Deallocate arrays 
00512 !
00513       deallocate (max1)
00514       deallocate (min1)
00515 
00516       deallocate (face2_max)
00517       deallocate (face2_min)
00518 
00519       deallocate (face1_max)
00520       deallocate (face1_min)
00521 !
00522 !===> All done
00523 !
00524 #ifdef VERBOSE
00525       print 9980, trim(ch_id), index, cyclic
00526 
00527       call psmile_flushstd
00528 #endif /* VERBOSE */
00529 !
00530 !  Formats:
00531 !
00532 9990 format (1x, a, ': psmile_get_cyclic_dir_3d_dble index', i2)
00533 9980 format (1x, a, ': psmile_get_cyclic_dir_3d_dble: eof index, cyclic', i2, l6)
00534 
00535       end subroutine PSMILe_Get_cyclic_dir_3d_dble