psmile_neigh_tricu_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_tricu_3d_reg
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_tricu_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_tricu_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 !     Location of indices for neighbors_3d
00067 !
00068 !           13   14   15   16
00069 !
00070 !     ^      9   10   11   12
00071 !     |
00072 !     |      5    6    7    8
00073 !
00074 !     j      1    2    3    4
00075 !
00076 !       i ---->
00077 
00078       Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)
00079 
00080 !     Indices of neighbor locations (interpolation bases)
00081 
00082       Integer, Intent (Out)           :: ierror
00083 
00084 !     Returns the error code of PSMILE_Neigh_tricu_3d_reg;
00085 !             ierror = 0 : No error
00086 !             ierror > 0 : Severe error
00087 !
00088 ! !DEFINED PARAMETERS:
00089 !
00090 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00091 !               = 2 ** ndim_3d
00092 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00093 !               = 2 ** ndim_2d
00094 !
00095       Integer, parameter              :: n_corners_3d = 2 * 2 * 2 * 8
00096       Integer, parameter              :: n_corners_2d = 2 * 2 * 4
00097       Integer, parameter              :: n_corners_1d = 2 * 2
00098 !
00099 ! !LOCAL VARIABLES
00100 !
00101 !     ... for locations searched
00102 !
00103       Integer                         :: i, j, k, n
00104       Integer                         :: nbeg, nend, nprev
00105       Integer                         :: lenij, lenijk, n_corners
00106 !
00107 !     ... For locations controlled
00108 !
00109       Integer                         :: length (ndim_3d)
00110 !
00111       Integer                         :: nca (ndim_3d)
00112 !
00113 !     ... for error handling
00114 !
00115       Integer, parameter              :: nerrp = 2
00116       Integer                         :: ierrp (nerrp)
00117 !
00118 !
00119 ! !DESCRIPTION:
00120 !
00121 ! Subroutine "PSMILe_Neigh_tricu_3d_reg" searches the "n_corners"
00122 ! (8, 4 or 2) neighbours (3d-indices of the interpolation bases) for the
00123 ! selected interpolation type (tri-, bi-linear or linear, respectively)
00124 ! on the the method-grid (x_coords, y_coords, z_coords)
00125 ! for the subgrid coords sent by the destination process.
00126 !
00127 ! Subroutine "PSMILe_Neigh_tricu_3d_reg" is designed for (target) grids
00128 ! of type "PRISM_Reglonlatvrt".
00129 !
00130 ! !REVISION HISTORY:
00131 !
00132 !   Date      Programmer   Description
00133 ! ----------  ----------   -----------
00134 ! 03.07.21    H. Ritzdorf  created trili
00135 ! 03.07.21    R. Redler    modified for cubic interpolation
00136 !
00137 !EOP
00138 !----------------------------------------------------------------------
00139 !
00140 !  $Id: psmile_neigh_tricu_3d_reg.F90 2687 2010-10-28 15:15:52Z coquart $
00141 !  $Author: coquart $
00142 !
00143    Character(len=len_cvs_string), save :: mycvs = 
00144        '$Id: psmile_neigh_tricu_3d_reg.F90 2687 2010-10-28 15:15:52Z coquart $'
00145 !
00146 !----------------------------------------------------------------------
00147 !
00148       data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00149 !
00150 !----------------------------------------------------------------------
00151 !
00152 !  Initialization
00153 !
00154 #ifdef VERBOSE
00155       print 9990, trim(ch_id)
00156 
00157       call psmile_flushstd
00158 #endif /* VERBOSE */
00159 !
00160 !===> Control interpolation mode
00161 !
00162       if (interpolation_mode < 1 .or. interpolation_mode > ndim_3d) then
00163          ierrp (1) = interpolation_mode
00164          ierror = PRISM_Error_Internal
00165 
00166          call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_tricu_irreg2', &
00167                              ierrp, 1, __FILE__, __LINE__ )
00168          return
00169       endif
00170 !
00171 !===> Initialization
00172 !
00173 !  lenij  = Number of entries in first 2 directions
00174 !  lenijk = Number of entries in all   3 directions
00175 !
00176       ierror = 0
00177 !
00178       nprev = npreva
00179       n_corners = nca (interpolation_mode)
00180 !
00181       lenij  = nlocs (1) * nlocs (2)
00182       lenijk = lenij * nlocs (3)
00183 !
00184       nbeg = nprev + 1
00185       nend = nprev + lenijk
00186 !
00187       length (1:ndim_3d) = grid_valid_shape(2,1:ndim_3d) - &
00188                            grid_valid_shape(1,1:ndim_3d) + 1
00189 !
00190 #ifdef PRISM_ASSERTION
00191 !
00192 !===> Is the dimension of array "neighbors_3d" sufficient
00193 !
00194       if (num_neigh < n_corners) then
00195          print 9970, trim(ch_id), num_neigh, n_corners
00196          call psmile_assert (__FILE__, __LINE__, &
00197                              'Number of neighbors is insufficient')
00198       endif
00199 !
00200 !===> 
00201 !
00202       if (nloc < nprev + nlocs (1) * nlocs (2) * nlocs (3)) then
00203          print *, trim(ch_id), " nloc, nprev, nlocs ", nloc, nprev, nlocs
00204          call psmile_assert (__FILE__, __LINE__, &
00205                              'nloc < nprev + PRODUCT (nlocs) ')
00206       endif
00207 !
00208 !===> Are the locations within the correct shape ?
00209 !     Note: the scrloc's are source locations in the 
00210 !           method grid which has an virtual cell 0.
00211 !
00212          do j = 1, ndim_3d
00213 !cdir vector
00214             do i = 1, nlocs (j)
00215                if ( srclocs(j)%vector(i) < grid_valid_shape(1,j) .or. &
00216                     srclocs(j)%vector(i) > grid_valid_shape(2,j)) exit
00217             end do
00218             if (i <= nlocs(j)) then
00219                print *, "Incorrect index in j, i", j, i, srclocs(j)%vector(i)
00220                call psmile_assert (__FILE__, __LINE__, 'wrong index')
00221             endif
00222          end do
00223 #endif
00224 !
00225 !===> Control special cases (2d hyperplanes)
00226 !
00227          do i = 1, ndim_3d
00228          if (length(i) == 1) then
00229                neighbors_3d (i, nbeg:nend, 1:n_corners) = grid_valid_shape(1,i)
00230             end if
00231          end do
00232 !
00233 !===> ... I indices of neighbour indices
00234 !
00235       if (length(1) > 1) then
00236 
00237             do j = 1, nlocs (2)
00238                neighbors_3d (1, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 6) = &
00239                     srclocs(1)%vector(:)
00240             end do
00241 
00242             do k = 2, nlocs (3)
00243                neighbors_3d (1, nprev+(k-1)*lenij+1:nprev+k*lenij, 6) = &
00244                     neighbors_3d (1, nprev+            1:nprev+  lenij, 6)
00245             end do
00246 
00247             do n = 0, 3
00248                neighbors_3d (1, nbeg:nend, 1+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) - 1
00249                neighbors_3d (1, nbeg:nend, 2+(n*4) ) = neighbors_3d (1, nbeg:nend, 6)
00250                neighbors_3d (1, nbeg:nend, 3+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) + 1
00251                neighbors_3d (1, nbeg:nend, 4+(n*4) ) = neighbors_3d (1, nbeg:nend, 6) + 2
00252             enddo
00253 
00254             if (n_corners > 16) &
00255                  neighbors_3d (1, nbeg:nend, 17:32) = neighbors_3d (1, nbeg:nend,  1:16)
00256 
00257          endif
00258 !
00259 !===> ... J indices of neighbour indices
00260 !
00261       if (length(2) > 1) then
00262 
00263             do j = 1, nlocs (2)
00264                neighbors_3d (2, nprev+(j-1)*nlocs(1)+1:nprev+j*nlocs(1), 6) = &
00265                     srclocs(2)%vector(j)
00266             end do
00267 
00268             do k = 2, nlocs (3)
00269                neighbors_3d (2, nprev+(k-1)*lenij+1:nprev+k*lenij, 6) = &
00270                neighbors_3d (2, nprev+            1:nprev+  lenij, 6)
00271             end do
00272 
00273             do n = 1, 4
00274                neighbors_3d (2, nbeg:nend,    n) = neighbors_3d (2, nbeg:nend, 6) - 1
00275                neighbors_3d (2, nbeg:nend,  4+n) = neighbors_3d (2, nbeg:nend, 6)
00276                neighbors_3d (2, nbeg:nend,  8+n) = neighbors_3d (2, nbeg:nend, 6) + 1
00277                neighbors_3d (2, nbeg:nend, 12+n) = neighbors_3d (2, nbeg:nend, 6) + 2
00278             end do
00279 
00280             if (n_corners > 16) then
00281                neighbors_3d (2, nbeg:nend, 17:20) = neighbors_3d (2, nbeg:nend,  1: 4)
00282                neighbors_3d (2, nbeg:nend, 21:24) = neighbors_3d (2, nbeg:nend,  5: 8)
00283                neighbors_3d (2, nbeg:nend, 25:28) = neighbors_3d (2, nbeg:nend,  9:12)
00284                neighbors_3d (2, nbeg:nend, 29:32) = neighbors_3d (2, nbeg:nend, 13:16)
00285             endif
00286 
00287          endif
00288 !
00289 !===> ... K indices of neighbour indices
00290 !
00291       if (length(3) > 1) then
00292 
00293             do k = 1, nlocs (3)
00294                neighbors_3d (3, nprev+(k-1)*lenij+1:nprev+k*lenij, 1) = &
00295                     srclocs(3)%vector(k)
00296             end do
00297 
00298             do n = 2, 16
00299                neighbors_3d (3, nbeg:nend, n) = neighbors_3d (3, nbeg:nend, 1)
00300             enddo
00301 
00302             if (n_corners > 16) then
00303                do n = 17, 32
00304                   neighbors_3d (3, nbeg:nend,n) = neighbors_3d (3, nbeg:nend, 1) + 1
00305                enddo
00306             endif
00307 
00308          endif
00309 !
00310 #if defined ( PRISM_ASSERTION ) && defined (TRANSF_SUP_2D_INTERP)
00311 !
00312 !===> Internal control
00313 !
00314          ic = nprev
00315 !
00316 kontrol: do k = 1, nlocs (3)
00317             do j  = 1, nlocs (2)
00318                do i = 1, nlocs (1)
00319                   ic = ic + 1
00320 
00321                   if ( neighbors_3d (1, ic, 6) /= srclocs(1)%vector(i)     .or. &
00322                        neighbors_3d (2, ic, 6) /= srclocs(2)%vector(j)     .or. &
00323                        neighbors_3d (3, ic, 6) /= srclocs(3)%vector(k)     .or. &
00324                        
00325                        neighbors_3d (1, ic, 7) /= srclocs(1)%vector(i) + 1 .or. &
00326                        neighbors_3d (2, ic, 7) /= srclocs(2)%vector(j)     .or. &
00327                        neighbors_3d (3, ic, 7) /= srclocs(3)%vector(k)     .or. &
00328                        
00329                        neighbors_3d (1, ic,11) /= srclocs(1)%vector(i) + 1 .or. &
00330                        neighbors_3d (2, ic,11) /= srclocs(2)%vector(j) + 1 .or. &
00331                        neighbors_3d (3, ic,11) /= srclocs(3)%vector(k)     .or. &
00332                        
00333                        neighbors_3d (1, ic,10) /= srclocs(1)%vector(i)     .or. &
00334                        neighbors_3d (2, ic,10) /= srclocs(2)%vector(j) + 1 .or. &
00335                        neighbors_3d (3, ic,10) /= srclocs(3)%vector(k) ) exit kontrol
00336 
00337                   if ( n_corners > 16 ) then
00338                      if ( neighbors_3d (1, ic,22) /= srclocs(1)%vector(i)     .or. &
00339                           neighbors_3d (2, ic,22) /= srclocs(2)%vector(j)     .or. &
00340                           neighbors_3d (3, ic,22) /= srclocs(3)%vector(k) + 1 .or. &
00341                           
00342                           neighbors_3d (1, ic,23) /= srclocs(1)%vector(i) + 1 .or. &
00343                           neighbors_3d (2, ic,23) /= srclocs(2)%vector(j)     .or. &
00344                           neighbors_3d (3, ic,23) /= srclocs(3)%vector(k) + 1 .or. &
00345                           
00346                           neighbors_3d (1, ic,27) /= srclocs(1)%vector(i)     .or. &
00347                           neighbors_3d (2, ic,27) /= srclocs(2)%vector(j) + 1 .or. &
00348                           neighbors_3d (3, ic,27) /= srclocs(3)%vector(k) + 1 .or. &
00349                           
00350                           neighbors_3d (1, ic,26) /= srclocs(1)%vector(i) + 1 .or. &
00351                           neighbors_3d (2, ic,26) /= srclocs(2)%vector(j) + 1 .or. &
00352                           neighbors_3d (3, ic,26) /= srclocs(3)%vector(k) + 1) exit kontrol
00353 
00354                   endif
00355 
00356                end do
00357             end do
00358          end do kontrol
00359 !
00360          if (i <= nlocs (1)) then
00361             print *, "Incorrect neighbour generated", i, j, k
00362             call psmile_assert (__FILE__, __LINE__, &
00363                                 'incorrect 3d neighbour generated')
00364          endif
00365 #endif /* PRISM_ASSERTION */
00366 !
00367 !===> ... set indices for cyclic coordinate directions
00368 !
00369          do j = 1, ndim_3d
00370          if (cyclic(j) .and. length(j) > 1) then
00371 !
00372                do n = 1, n_corners
00373 !cdir vector
00374                   do i = nbeg, nend
00375                      
00376                      if  ( neighbors_3d (j, i, n) > grid_valid_shape(2,j) ) then
00377 
00378                            neighbors_3d (j, i, n) = &
00379                            neighbors_3d (j, i, n) - length(j)
00380                      else if ( neighbors_3d (j, i, n) < grid_valid_shape(1,j) ) then
00381 
00382                            neighbors_3d (j, i, n) = &
00383                            neighbors_3d (j, i, n) + length(j)
00384                      endif
00385 
00386                   end do
00387 
00388                end do ! n
00389             endif
00390          end do ! j
00391 !
00392 !===> All done
00393 !
00394 #ifdef VERBOSE
00395       print 9980, trim(ch_id)
00396 
00397       call psmile_flushstd
00398 #endif /* VERBOSE */
00399 !
00400 !  Formats:
00401 !
00402 9990 format (1x, a, ': psmile_neigh_tricu_3d_reg')
00403 9980 format (1x, a, ': psmile_neigh_tricu_3d_reg: eof')
00404 
00405 #ifdef PRISM_ASSERTION
00406 9970 format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
00407 #endif
00408 
00409       end subroutine psmile_neigh_tricu_3d_reg