psmile_neigh_trili_3d_reg.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_Neigh_trili_3d_reg
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_trili_3d_reg (                 &
00012                grid_valid_shape, interpolation_mode, cyclic, &
00013                srclocs, nlocs, nloc, npreva, neighbors_3d,   &
00014                num_neigh, ierror)
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019 !
00020       use PSMILe, dummy_interface => PSMILe_Neigh_trili_3d_reg
00021 
00022       implicit none
00023 !
00024 ! !INPUT PARAMETERS:
00025 !
00026       Integer, Intent (In)            :: nloc
00027 !
00028 !     Total number of locations to be transferred
00029 !
00030       Integer, Intent (In)            :: nlocs (ndim_3d)
00031 !
00032 !     Number of locations in each direction
00033 !
00034       Type (integer_vector), Intent (In) :: srclocs (ndim_3d)
00035 !
00036 !     Indices of the (method) grid cell in which the point was found.
00037 
00038       Integer, Intent (In)            :: npreva
00039 !
00040 !     Number of locations for which neighbours were already searched.
00041 !
00042       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00043 
00044 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00045 !
00046       Integer,            Intent (In) :: interpolation_mode
00047 !
00048 !     Specifies the interpolation mode with
00049 !        interpolation_mode = 3 : Tri-linear or bilinear/linear interpolation
00050 !                                 (determines 8 corners)
00051 !        interpolation_mode = 2 : Bilinear interpolation
00052 !                                 (determines 4 corners)
00053 !        interpolation_mode = 1 : Linear   interpolation
00054 !                                 (determines 2 corners)
00055 !
00056       Logical,            Intent (In) :: cyclic (ndim_3d)
00057 !    
00058 !     Does the block have cyclic coordinates in the corresponding direction ?
00059 !
00060       Integer,            Intent (In) :: num_neigh
00061 !
00062 !     Last dimension of neighbors array "neighbors_3d"
00063 !
00064 ! !OUTPUT PARAMETERS:
00065 !
00066       Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)
00067 
00068 !     Indices of neighbor locations (interpolation bases)
00069 
00070       Integer, Intent (Out)           :: ierror
00071 
00072 !     Returns the error code of PSMILE_Neigh_trili_3d_reg;
00073 !             ierror = 0 : No error
00074 !             ierror > 0 : Severe error
00075 !
00076 ! !DEFINED PARAMETERS:
00077 !
00078 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00079 !               = 2 ** ndim_3d
00080 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00081 !               = 2 ** ndim_2d
00082 !
00083       Integer, Parameter              :: n_corners_3d = 2 * 2 * 2
00084       Integer, Parameter              :: n_corners_2d = 2 * 2
00085       Integer, Parameter              :: n_corners_1d = 2
00086 !
00087 ! !LOCAL VARIABLES
00088 !
00089 !     ... for locations searched
00090 !
00091 !     These long character strings cause problems for some preprocessing tools.
00092 !     as they cannot be properly resolved.
00093 !
00094 #if defined ( PRISM_ASSERTION ) && defined (TRANSF_SUP_2D_INTERP_HUHU)
00095       Integer                         :: ic
00096 #endif
00097       Integer                         :: i, j, k, n, nbeg, nend, nprev
00098       Integer                         :: lenij, lenijk, n_corners
00099 !
00100 !     ... For locations controlled
00101 !
00102       Integer                         :: length (ndim_3d)
00103 !
00104       Integer                         :: nca (ndim_3d)
00105 !
00106 !     ... for error handling
00107 !
00108       Integer, Parameter              :: nerrp = 2
00109       Integer                         :: ierrp (nerrp)
00110 !
00111 !
00112 ! !DESCRIPTION:
00113 !
00114 ! Subroutine "PSMILe_Neigh_trili_3d_reg" searches the "n_corners"
00115 ! (8, 4 or 2) neighbours (3d-indices of the interpolation bases) for the
00116 ! selected interpolation type (tri-, bi-linear or linear, respectively)
00117 ! on the the method-grid (x_coords, y_coords, z_coords)
00118 ! for the subgrid coords sent by the destination process.
00119 !
00120 ! Subroutine "PSMILe_Neigh_trili_3d_reg" is designed for (target) grids
00121 ! of type "PRISM_Reglonlatvrt".
00122 !
00123 ! !REVISION HISTORY:
00124 !
00125 !   Date      Programmer   Description
00126 ! ----------  ----------   -----------
00127 ! 03.07.05    H. Ritzdorf  created
00128 !
00129 !EOP
00130 !----------------------------------------------------------------------
00131 !
00132 !  $Id: psmile_neigh_trili_3d_reg.F90 2687 2010-10-28 15:15:52Z coquart $
00133 !  $Author: coquart $
00134 !
00135    Character(len=len_cvs_string), save :: mycvs = 
00136        '$Id: psmile_neigh_trili_3d_reg.F90 2687 2010-10-28 15:15:52Z coquart $'
00137 !
00138 !----------------------------------------------------------------------
00139 !
00140       data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00141 !
00142 !----------------------------------------------------------------------
00143 !
00144 !  Initialization
00145 !
00146 #ifdef VERBOSE
00147       print 9990, trim(ch_id)
00148 
00149       call psmile_flushstd
00150 #endif /* VERBOSE */
00151 !
00152 !===> Control interpolation mode
00153 !
00154       if (interpolation_mode < 1 .or. interpolation_mode > ndim_3d) then
00155          ierrp (1) = interpolation_mode
00156          ierror = PRISM_Error_Internal
00157 
00158          call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_trili_reg', &
00159                              ierrp, 1, __FILE__, __LINE__ )
00160          return
00161       endif
00162 !
00163 !===> Initialization
00164 !
00165 !  lenij  = Number of entries in first 2 directions
00166 !  lenijk = Number of entries in all   3 directions
00167 !
00168       ierror = 0
00169 !
00170       nprev = npreva
00171       n_corners = nca (interpolation_mode)
00172 !
00173       lenij  = nlocs (1) * nlocs (2)
00174       lenijk = lenij * nlocs (3)
00175 !
00176       nbeg = nprev + 1
00177       nend = nprev + lenijk
00178 !
00179       length (1:ndim_3d) = grid_valid_shape(2,1:ndim_3d) - &
00180                            grid_valid_shape(1,1:ndim_3d) + 1
00181 !
00182 #ifdef PRISM_ASSERTION
00183 !
00184 !===> Is the dimension of array "neighbors_3d" sufficient
00185 !
00186       if (num_neigh < n_corners) then
00187          print 9970, trim(ch_id), num_neigh, n_corners
00188          call psmile_assert (__FILE__, __LINE__, &
00189                              'Number of neighbors is insufficient')
00190       endif
00191 !
00192 !===> 
00193 !
00194       if (nloc < nprev + nlocs (1) * nlocs (2) * nlocs (3)) then
00195          print *, trim(ch_id), " nloc, nprev, nlocs ", nloc, nprev, nlocs
00196          call psmile_assert (__FILE__, __LINE__, &
00197                              'nloc < nprev + PRODUCT (nlocs) ')
00198       endif
00199 !
00200 !===> Are the locations within the correct shape ?
00201 !     Note: the scrloc's are source locations in the 
00202 !           method grid which has an virtual cell 0.
00203 !
00204          do j = 1, ndim_3d
00205 !cdir vector
00206             do i = 1, nlocs (j)
00207 !
00208             if (srclocs(j)%vector(i) < grid_valid_shape(1,j)-1 .or. &
00209                 srclocs(j)%vector(i) > grid_valid_shape(2,j))  exit
00210             end do
00211 !
00212          if (i <= nlocs(j)) then
00213             print *, "Incorrect index in j, i", j, i, srclocs(j)%vector(i)
00214             call psmile_assert (__FILE__, __LINE__, &
00215                                 'wrong index')
00216          endif
00217          end do
00218 #endif
00219 !
00220 !===> Generate entries (assuming mask is not available or all
00221 !     mask points are defined)
00222 !
00223 !     ic = nprev
00224 !        do k = 1, nlocs (3)
00225 !           do j  = 1, nlocs (2)
00226 !              do i = 1, nlocs (1)
00227 !              ic = ic + 1
00228 ! 
00229 !              neighbors_3d (1, ic, 1) = srclocs(1)%vector(i)
00230 !              neighbors_3d (2, ic, 1) = srclocs(2)%vector(j)
00231 !              neighbors_3d (3, ic, 1) = srclocs(3)%vector(k)
00232 !
00233 !              neighbors_3d (1, ic, 2) = srclocs(1)%vector(i) + 1
00234 !              neighbors_3d (2, ic, 2) = srclocs(2)%vector(j)
00235 !              neighbors_3d (3, ic, 2) = srclocs(3)%vector(k)
00236 !
00237 ! expected by the transformer
00238 !              neighbors_3d (1, ic, 3) = srclocs(1)%vector(i) + 1
00239 !              neighbors_3d (2, ic, 3) = srclocs(2)%vector(j) + 1
00240 !              neighbors_3d (3, ic, 3) = srclocs(3)%vector(k)
00241 !
00242 !              neighbors_3d (1, ic, 4) = srclocs(1)%vector(i)
00243 !              neighbors_3d (2, ic, 4) = srclocs(2)%vector(j) + 1
00244 !              neighbors_3d (3, ic, 4) = srclocs(3)%vector(k)
00245 !
00246 !              neighbors_3d (1, ic, 5) = srclocs(1)%vector(i)
00247 !              neighbors_3d (2, ic, 5) = srclocs(2)%vector(j)
00248 !              neighbors_3d (3, ic, 5) = srclocs(3)%vector(k) + 1
00249 !
00250 !              neighbors_3d (1, ic, 6) = srclocs(1)%vector(i) + 1
00251 !              neighbors_3d (2, ic, 6) = srclocs(2)%vector(j)
00252 !              neighbors_3d (3, ic, 6) = srclocs(3)%vector(k) + 1
00253 !
00254 ! expected by the transformer
00255 !              neighbors_3d (1, ic, 7) = srclocs(1)%vector(i) + 1
00256 !              neighbors_3d (2, ic, 7) = srclocs(2)%vector(j) + 1
00257 !              neighbors_3d (3, ic, 7) = srclocs(3)%vector(k) + 1
00258 !
00259 !              neighbors_3d (1, ic, 8) = srclocs(1)%vector(i)
00260 !              neighbors_3d (2, ic, 8) = srclocs(2)%vector(j) + 1
00261 !              neighbors_3d (3, ic, 8) = srclocs(3)%vector(k) + 1
00262 !              end do
00263 !           end do
00264 !        end do
00265 !
00266 !===> Control special cases (2d hyperplanes)
00267 !
00268          do i = 1, ndim_3d
00269          if (length(i) == 1) then
00270             neighbors_3d (i, nbeg:nend, 1:n_corners) = grid_valid_shape(1,i)
00271          endif
00272          end do
00273 !
00274 !===> ... I indices of neighbour indices
00275 !
00276       if (length(1) > 1) then
00277 !
00278 
00279 #   define TRANSF_SUP_2D_INTERP
00280 !
00281 #ifdef TRANSF_SUP_2D_INTERP
00282 !
00283 ! Note: The transformer doesn't currently perform 2d interpolation
00284 !       for a point which lies on the face of the "last cell"
00285 
00286             do j = 1, nlocs (2)
00287             neighbors_3d (1, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00288                srclocs(1)%vector(:)
00289             end do
00290 #else
00291          if ( cyclic (1) ) then
00292                do j = 1, nlocs (2)
00293                neighbors_3d (1, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00294                   srclocs(1)%vector(:)
00295                end do
00296          else
00297 
00298                do j = 1, nlocs (2)
00299                neighbors_3d (1, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00300                   min (srclocs(1)%vector(:), grid_valid_shape(2,1)-1)
00301                end do
00302          endif
00303 #endif
00304 !
00305             do k = 2, nlocs (3)
00306             neighbors_3d (1, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00307             neighbors_3d (1, nprev+            1:nprev+  lenij, 1)
00308             end do
00309 !
00310          neighbors_3d (1, nbeg:nend, 2) = neighbors_3d (1, nbeg:nend, 1) + 1
00311 !
00312          if (n_corners > 2) then
00313             neighbors_3d (1, nbeg:nend, 3) = neighbors_3d (1, nbeg:nend, 2)
00314             neighbors_3d (1, nbeg:nend, 4) = neighbors_3d (1, nbeg:nend, 1)
00315 
00316             if (n_corners > 4) then
00317                neighbors_3d (1, nbeg:nend, 5) = neighbors_3d (1, nbeg:nend, 1)
00318                neighbors_3d (1, nbeg:nend, 8) = neighbors_3d (1, nbeg:nend, 1)
00319                neighbors_3d (1, nbeg:nend, 6) = neighbors_3d (1, nbeg:nend, 2)
00320                neighbors_3d (1, nbeg:nend, 7) = neighbors_3d (1, nbeg:nend, 2)
00321             endif
00322          endif
00323       endif
00324 !
00325 !===> ... J indices of neighbour indices
00326 !
00327       if (length(2) > 1) then
00328 !
00329 #ifdef TRANSF_SUP_2D_INTERP
00330             do j = 1, nlocs (2)
00331             neighbors_3d (2, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00332                srclocs(2)%vector(j)
00333             end do
00334 #else
00335          if ( cyclic (2) ) then
00336 
00337             do j = 1, nlocs (2)
00338             neighbors_3d (2, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00339                srclocs(2)%vector(j)
00340             end do
00341 
00342          else
00343 
00344             do j = 1, nlocs (2)
00345             neighbors_3d (2, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 1) = &
00346                min(srclocs(2)%vector(j), grid_valid_shape(2,2)-1)
00347             end do
00348 
00349          endif
00350 #endif
00351 !
00352             do k = 2, nlocs (3)
00353             neighbors_3d (2, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00354             neighbors_3d (2, nprev+            1:nprev+  lenij, 1)
00355             end do
00356 !
00357          neighbors_3d (2, nbeg:nend, 2) = neighbors_3d (2, nbeg:nend, 1)
00358 !
00359          if (n_corners > 2) then
00360             neighbors_3d (2, nbeg:nend, 3) = neighbors_3d (2, nbeg:nend, 1) + 1
00361             neighbors_3d (2, nbeg:nend, 4) = neighbors_3d (2, nbeg:nend, 3)
00362 
00363             if (n_corners > 4) then
00364                neighbors_3d (2, nbeg:nend, 5) = neighbors_3d (2, nbeg:nend, 1)
00365                neighbors_3d (2, nbeg:nend, 6) = neighbors_3d (2, nbeg:nend, 1)
00366 !
00367                neighbors_3d (2, nbeg:nend, 7) = neighbors_3d (2, nbeg:nend, 3)
00368                neighbors_3d (2, nbeg:nend, 8) = neighbors_3d (2, nbeg:nend, 3)
00369             endif
00370          endif
00371       endif
00372 !
00373 !===> ... K indices of neighbour indices
00374 !
00375       if (length(3) > 1) then
00376 
00377 #ifdef TRANSF_SUP_2D_INTERP
00378             do k = 1, nlocs (3)
00379             neighbors_3d (3, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00380                srclocs(3)%vector(k)
00381             end do
00382 #else
00383          if ( cyclic (3) .or. n_corners == n_corners_2d) then
00384             do k = 1, nlocs (3)
00385             neighbors_3d (3, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00386                srclocs(3)%vector(k)
00387             end do
00388 
00389          else
00390 
00391             do k = 1, nlocs (3)
00392             neighbors_3d (3, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00393                min (srclocs(3)%vector(k), grid_valid_shape(2, 3) - 1)
00394             end do
00395 
00396          endif
00397 #endif
00398 !
00399          neighbors_3d (3, nbeg:nend, 2) = neighbors_3d (3, nbeg:nend, 1)
00400 !
00401          if (n_corners > 2) then
00402             neighbors_3d (3, nbeg:nend, 3) = neighbors_3d (3, nbeg:nend, 1)
00403             neighbors_3d (3, nbeg:nend, 4) = neighbors_3d (3, nbeg:nend, 1)
00404 !
00405             if (n_corners > 4) then
00406                neighbors_3d (3, nbeg:nend, 5) = neighbors_3d (3, nbeg:nend, 1) &
00407                                                 + 1
00408                neighbors_3d (3, nbeg:nend, 6) = neighbors_3d (3, nbeg:nend, 5)
00409                neighbors_3d (3, nbeg:nend, 7) = neighbors_3d (3, nbeg:nend, 5)
00410                neighbors_3d (3, nbeg:nend, 8) = neighbors_3d (3, nbeg:nend, 5)
00411             endif
00412          endif
00413 
00414       endif
00415 !
00416 #if defined ( PRISM_ASSERTION ) && defined (TRANSF_SUP_2D_INTERP_HUHU)
00417 !
00418 !===> Internal control
00419 !
00420       ic = nprev
00421 !
00422 kontrol: do k = 1, nlocs (3)
00423             do j  = 1, nlocs (2)
00424                do i = 1, nlocs (1)
00425                ic = ic + 1
00426  
00427                if (neighbors_3d (1, ic, 1) /= srclocs(1)%vector(i)     .or. &
00428                    neighbors_3d (2, ic, 1) /= srclocs(2)%vector(j)     .or. &
00429                    neighbors_3d (3, ic, 1) /= srclocs(3)%vector(k)     .or. &
00430                    neighbors_3d (1, ic, 2) /= srclocs(1)%vector(i) + 1 .or. &
00431                    neighbors_3d (2, ic, 2) /= srclocs(2)%vector(j)     .or. &
00432                    neighbors_3d (3, ic, 2) /= srclocs(3)%vector(k)     .or. &
00433                    neighbors_3d (1, ic, 3) /= srclocs(1)%vector(i) + 1 .or. &
00434                    neighbors_3d (2, ic, 3) /= srclocs(2)%vector(j) + 1 .or. &
00435                    neighbors_3d (3, ic, 3) /= srclocs(3)%vector(k)     .or. &
00436                    neighbors_3d (1, ic, 4) /= srclocs(1)%vector(i)     .or. &
00437                    neighbors_3d (2, ic, 4) /= srclocs(2)%vector(j) + 1 .or. &
00438                    neighbors_3d (3, ic, 4) /= srclocs(3)%vector(k)     .or. &
00439                    neighbors_3d (1, ic, 5) /= srclocs(1)%vector(i)     .or. &
00440                    neighbors_3d (2, ic, 5) /= srclocs(2)%vector(j)     .or. &
00441                    neighbors_3d (3, ic, 5) /= srclocs(3)%vector(k) + 1 .or. &
00442                    neighbors_3d (1, ic, 6) /= srclocs(1)%vector(i) + 1 .or. &
00443                    neighbors_3d (2, ic, 6) /= srclocs(2)%vector(j)     .or. &
00444                    neighbors_3d (3, ic, 6) /= srclocs(3)%vector(k) + 1 .or. &
00445                    neighbors_3d (1, ic, 7) /= srclocs(1)%vector(i) + 1 .or. &
00446                    neighbors_3d (2, ic, 7) /= srclocs(2)%vector(j) + 1 .or. &
00447                    neighbors_3d (3, ic, 7) /= srclocs(3)%vector(k) + 1 .or. &
00448                    neighbors_3d (1, ic, 8) /= srclocs(1)%vector(i)     .or. &
00449                    neighbors_3d (2, ic, 8) /= srclocs(2)%vector(j) + 1 .or. &
00450                    neighbors_3d (3, ic, 8) /= srclocs(3)%vector(k) + 1) exit kontrol
00451                end do
00452             end do
00453          end do kontrol
00454 !
00455          if (i <= nlocs (1)) then
00456             print *, "Incorrect neighbour generated", i, j, k
00457             call psmile_assert (__FILE__, __LINE__, &
00458                                 'incorrect 3d neighbour generated')
00459          endif
00460 #endif /* PRISM_ASSERTION */
00461 !
00462 !===> ... set indices for cyclic coordinate directions
00463 !
00464             do j = 1, ndim_3d
00465             if (cyclic(j) .and. length(j) > 1) then
00466 !
00467                   do n = 1, n_corners
00468 !cdir vector
00469                      do i = nbeg, nend
00470                      if (neighbors_3d (j, i, n) < grid_valid_shape(1,j)) then
00471 
00472                          neighbors_3d (j, i, n) = &
00473                          neighbors_3d (j, i, n) + length (j)
00474 
00475                      else if ( neighbors_3d (j, i, n) > grid_valid_shape(2,j)) then
00476                          neighbors_3d (j, i, n) = &
00477                          neighbors_3d (j, i, n) - length (j)
00478 
00479                      endif
00480                      end do
00481                   end do ! n
00482             endif
00483             end do ! j
00484 !
00485 !===> All done
00486 !
00487 #ifdef VERBOSE
00488       print 9980, trim(ch_id)
00489 
00490       call psmile_flushstd
00491 #endif /* VERBOSE */
00492 !
00493 !  Formats:
00494 !
00495 9990 format (1x, a, ': psmile_neigh_trili_3d_reg')
00496 9980 format (1x, a, ': psmile_neigh_trili_3d_reg: eof')
00497 
00498 #ifdef PRISM_ASSERTION
00499 9970 format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
00500 #endif
00501 
00502       end subroutine psmile_neigh_trili_3d_reg