psmile_neigh_trili_irreg2.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_irreg2
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_trili_irreg2 (                 &
00012                grid_valid_shape, interpolation_mode, cyclic, &
00013                srcloc, srclocz, nlocs, nloc, npreva,         &
00014                neighbors_3d, num_neigh, ierror)
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019 !
00020       use PSMILe, dummy_interface => PSMILe_Neigh_trili_irreg2
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_2d)
00031 !
00032 !     Number of locations in each direction
00033 !
00034       Integer, Intent (In)            :: srcloc  (ndim_2d, nlocs (1))
00035       Integer, Intent (In)            :: srclocz (nlocs (2))
00036 !
00037 !     Indices of the (method) grid cell in which the point was found.
00038 
00039       Integer, Intent (In)            :: npreva
00040 !
00041 !     Number of locations for which neighbours were already searched.
00042 !
00043       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00044 
00045 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00046 !
00047       Integer,            Intent (In) :: interpolation_mode
00048 !
00049 !     Specifies the interpolation mode with
00050 !        interpolation_mode = 3 : Tri-linear or bilinear/linear interpolation
00051 !                                 (determines 8 corners)
00052 !        interpolation_mode = 2 : Bilinear interpolation
00053 !                                 (determines 4 corners)
00054 !        interpolation_mode = 1 : Linear   interpolation
00055 !                                 (determines 2 corners)
00056 !
00057       Logical,            Intent (In) :: cyclic (ndim_3d)
00058 !    
00059 !     Does the block have cyclic coordinates in the corresponding direction ?
00060 !
00061       Integer,            Intent (In) :: num_neigh
00062 !
00063 !     Last dimension of neighbors array "neighbors_3d"
00064 !
00065 ! !OUTPUT PARAMETERS:
00066 !
00067       Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)
00068 
00069 !     Indices of neighbor locations (interpolation bases)
00070 
00071       Integer, Intent (Out)           :: ierror
00072 
00073 !     Returns the error code of PSMILE_Neigh_trili_irreg2;
00074 !             ierror = 0 : No error
00075 !             ierror > 0 : Severe error
00076 !
00077 ! !DEFINED PARAMETERS:
00078 !
00079 !  N_CORNERS_3D = Number of corners for a 3-d interpolation
00080 !               = 2 ** ndim_3d
00081 !  N_CORNERS_2D = Number of corners for a 2-d interpolation
00082 !               = 2 ** ndim_2d
00083 !
00084       Integer, parameter              :: n_corners_3d = 2 * 2 * 2
00085       Integer, parameter              :: n_corners_2d = 2 * 2
00086       Integer, parameter              :: n_corners_1d = 2
00087 !
00088 ! !LOCAL VARIABLES
00089 !
00090 !     ... for locations searched
00091 !
00092       Integer                         :: i, j, k, n
00093       Integer                         :: nbeg, nend, nprev
00094       Integer                         :: lenijk, n_corners
00095 !
00096 !     ... For locations controlled
00097 !
00098       Integer                         :: length (ndim_3d)
00099 !
00100       Integer                         :: nca (ndim_3d)
00101 !
00102 !     ... for error handling
00103 !
00104       Integer, parameter              :: nerrp = 2
00105       Integer                         :: ierrp (nerrp)
00106 !
00107 !
00108 ! !DESCRIPTION:
00109 !
00110 ! Subroutine "PSMILe_Neigh_trili_irreg2" searches the "n_corners"
00111 ! (8, 4 or 2) neighbours (3d-indices of the interpolation bases) for the
00112 ! selected interpolation type (tri-, bi-linear or linear, respectively)
00113 ! on the the method-grid (x_coords, y_coords, z_coords)
00114 ! for the subgrid coords sent by the destination process.
00115 !
00116 ! Subroutine "PSMILe_Neigh_trili_irreg2" is designed for (target) grids
00117 ! of type "PRISM_Irrlonlat_regvrt".
00118 !
00119 ! !REVISION HISTORY:
00120 !
00121 !   Date      Programmer   Description
00122 ! ----------  ----------   -----------
00123 ! 03.07.21    H. Ritzdorf  created
00124 !
00125 !EOP
00126 !----------------------------------------------------------------------
00127 !
00128 !  $Id: psmile_neigh_trili_irreg2.F90 2687 2010-10-28 15:15:52Z coquart $
00129 !  $Author: coquart $
00130 !
00131    Character(len=len_cvs_string), save :: mycvs = 
00132        '$Id: psmile_neigh_trili_irreg2.F90 2687 2010-10-28 15:15:52Z coquart $'
00133 !
00134 !----------------------------------------------------------------------
00135 !
00136       data nca /n_corners_1d, n_corners_2d, n_corners_3d/
00137 !
00138 !----------------------------------------------------------------------
00139 !
00140 !  Initialization
00141 !
00142 #ifdef VERBOSE
00143       print 9990, trim(ch_id)
00144 
00145       call psmile_flushstd
00146 #endif /* VERBOSE */
00147 !
00148 !===> Control interpolation mode
00149 !
00150       if (interpolation_mode < 1 .or. interpolation_mode > ndim_3d) then
00151          ierrp (1) = interpolation_mode
00152          ierror = PRISM_Error_Internal
00153 
00154          call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_trili_irreg2', &
00155                              ierrp, 1, __FILE__, __LINE__ )
00156          return
00157       endif
00158 !
00159 !===> Initialization
00160 !
00161 !  lenijk = Number of entries in all   3 directions
00162 !
00163       ierror = 0
00164 !
00165       nprev = npreva
00166       n_corners = nca (interpolation_mode)
00167 !
00168       lenijk = nlocs (1) * nlocs (2)
00169       nbeg = nprev + 1
00170       nend = nprev + lenijk
00171 !
00172       length (1:ndim_3d) = grid_valid_shape(2,1:ndim_3d) - &
00173                            grid_valid_shape(1,1:ndim_3d) + 1
00174 !
00175 #ifdef PRISM_ASSERTION
00176 !
00177 !===> Is the dimension of array "neighbors_3d" sufficient
00178 !
00179       if (num_neigh < n_corners) then
00180          print 9970, trim(ch_id), num_neigh, n_corners
00181          call psmile_assert (__FILE__, __LINE__, &
00182                              'Number of neighbors is insufficient')
00183       endif
00184 !
00185 !===> 
00186 !
00187       if (nloc < nprev + nlocs (1) * nlocs (2) ) then
00188          print *, trim(ch_id), " nloc, nprev, nlocs ", nloc, nprev, nlocs
00189          call psmile_assert (__FILE__, __LINE__, &
00190                              'nloc < nprev + PRODUCT (nlocs) ')
00191       endif
00192 !
00193 !===> Are the locations within the correct shape ?
00194 !     Note: the scrloc's are source locations in the 
00195 !           method grid which has an virtual cell 0.
00196 !
00197 !cdir vector
00198          do i = 1, nlocs (1)
00199             if ( srcloc(1,i) < grid_valid_shape(1,1) -1 .or. &
00200                  srcloc(1,i) > grid_valid_shape(2,1)    .or. &
00201                  srcloc(2,i) < grid_valid_shape(1,2) -1 .or. &
00202                  srcloc(2,i) > grid_valid_shape(2,2)) exit
00203          end do
00204 
00205       if (i <= nlocs(1)) then
00206          print *, "Incorrect index in srcloc, i", i, srcloc(:, i)
00207          call psmile_assert (__FILE__, __LINE__, 'wrong index')
00208       endif
00209 
00210 !cdir vector
00211          do i = 1, nlocs (2)
00212             if ( srclocz(i) < grid_valid_shape(1,3)-1 .or. &
00213                  srclocz(i) > grid_valid_shape(2,3)) exit
00214          end do
00215 !
00216          if (i <= nlocs(2)) then
00217             print *, "Incorrect index in srclocs(2), i", i, srclocz(i), &
00218                      grid_valid_shape(1,3)-1, grid_valid_shape(2,3)
00219             call psmile_assert (__FILE__, __LINE__, 'wrong index')
00220          endif
00221 #endif
00222 !
00223 !===> Control special cases (2d hyperplanes)
00224 !
00225          do i = 1, ndim_3d
00226             if (length(i) == 1) then
00227                neighbors_3d (i, nbeg:nend, 1:n_corners) = grid_valid_shape(1,i)
00228             endif
00229          end do
00230 !
00231 !===> Generate entries (assuming mask is not available or all
00232 !     mask points are defined)
00233 !
00234 #define TRANSF_SUP_2D_INTERP
00235 
00236 #ifdef TRANSF_SUP_2D_INTERP
00237 !
00238 ! Note: The transformer doesn't currently perform 2d interpolation
00239 !       for a point which lies on the face of the "last cell"
00240 
00241       neighbors_3d (1:2, nprev+1:nprev+nlocs(1), 1) = srcloc (1:2, 1:nlocs(1))
00242 #else
00243       if (length(1) > 1) then
00244          if (cyclic (1)) then
00245             neighbors_3d (1, nprev+1:nprev+nlocs(1), 1) = srcloc (1, 1:nlocs(1))
00246 
00247          else
00248 
00249             neighbors_3d (1, nprev+1:nprev+nlocs(1), 1) = &
00250                  min (srcloc (1, 1:nlocs(1)), grid_valid_shape(2,1)-1)
00251          endif
00252       endif
00253 
00254       if (length(2) > 1) then
00255          if (cyclic (2)) then
00256             neighbors_3d (2, nprev+1:nprev+nlocs(1), 1) = srcloc (2, 1:nlocs(1))
00257 
00258          else
00259             neighbors_3d (2, nprev+1:nprev+nlocs(1), 1) = &
00260                  min (srcloc (2, 1:nlocs(1)), grid_valid_shape(2,2)-1)
00261          endif
00262       endif
00263 #endif
00264 !
00265 !===> Duplicate first 2 indices for third directions
00266 !
00267          do k = 2, nlocs (2)
00268          neighbors_3d (1:2, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 1) = &
00269          neighbors_3d (1:2, nprev+               1:nprev+  nlocs(1), 1)
00270          end do
00271 !
00272 !===> 
00273 !
00274 #ifdef TRANSF_SUP_2D_INTERP
00275          do k = 1, nlocs (2)
00276          neighbors_3d (3, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 1) = srclocz (k)     
00277          end do
00278 #else
00279       if (length(3) > 1) then
00280          if (cyclic (3) .or. n_corners == n_corners_2d) then
00281                do k = 1, nlocs (2)
00282                 neighbors_3d (3, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 1) = &
00283                        srclocz (k)     
00284                end do
00285 
00286           else
00287 
00288                do k = 1, nlocs (2)
00289                 neighbors_3d (3, nprev+(k-1)*nlocs(1)+1:nprev+k*nlocs(1), 1) = &
00290                        min (srclocz (k), grid_valid_shape(2,3)-1)
00291                end do
00292 
00293           endif
00294        endif
00295 #endif
00296 !
00297 !===> ... I indices of neighbour indices
00298 !
00299        if (length(1) > 1) then
00300           neighbors_3d (1, nbeg:nend, 2) = neighbors_3d (1, nbeg:nend, 1) + 1
00301 !
00302          if (n_corners > 2) then
00303             neighbors_3d (1, nbeg:nend, 3) = neighbors_3d (1, nbeg:nend, 2)
00304             neighbors_3d (1, nbeg:nend, 4) = neighbors_3d (1, nbeg:nend, 1)
00305 
00306             if (n_corners > 4) then
00307                neighbors_3d (1, nbeg:nend, 5) = neighbors_3d (1, nbeg:nend, 1)
00308                neighbors_3d (1, nbeg:nend, 8) = neighbors_3d (1, nbeg:nend, 1)
00309                neighbors_3d (1, nbeg:nend, 6) = neighbors_3d (1, nbeg:nend, 2)
00310                neighbors_3d (1, nbeg:nend, 7) = neighbors_3d (1, nbeg:nend, 2)
00311             endif
00312          endif
00313       endif
00314 !
00315 !===> ... J indices of neighbour indices
00316 !
00317       if (length(2) > 1) then
00318          neighbors_3d (2, nbeg:nend, 2) = neighbors_3d (2, nbeg:nend, 1)
00319 !
00320          if (n_corners > 2) then
00321             neighbors_3d (2, nbeg:nend, 3) = neighbors_3d (2, nbeg:nend, 1) + 1
00322             neighbors_3d (2, nbeg:nend, 4) = neighbors_3d (2, nbeg:nend, 3)
00323 
00324             if (n_corners > 4) then
00325                neighbors_3d (2, nbeg:nend, 5) = neighbors_3d (2, nbeg:nend, 1)
00326                neighbors_3d (2, nbeg:nend, 6) = neighbors_3d (2, nbeg:nend, 1)
00327 !
00328                neighbors_3d (2, nbeg:nend, 7) = neighbors_3d (2, nbeg:nend, 3)
00329                neighbors_3d (2, nbeg:nend, 8) = neighbors_3d (2, nbeg:nend, 3)
00330             endif
00331          endif
00332       endif
00333 !
00334 !===> ... K indices of neighbour indices
00335 !
00336       if (length(3) > 1) then
00337          neighbors_3d (3, nbeg:nend, 2) = neighbors_3d (3, nbeg:nend, 1)
00338 
00339          if (n_corners > 2) then
00340             neighbors_3d (3, nbeg:nend, 3) = neighbors_3d (3, nbeg:nend, 1)
00341             neighbors_3d (3, nbeg:nend, 4) = neighbors_3d (3, nbeg:nend, 1)
00342 
00343             if (n_corners > 4) then
00344                neighbors_3d (3, nbeg:nend, 5) = neighbors_3d (3, nbeg:nend, 1) + 1
00345                neighbors_3d (3, nbeg:nend, 6) = neighbors_3d (3, nbeg:nend, 5)
00346                neighbors_3d (3, nbeg:nend, 7) = neighbors_3d (3, nbeg:nend, 5)
00347                neighbors_3d (3, nbeg:nend, 8) = neighbors_3d (3, nbeg:nend, 5)
00348             endif
00349          endif
00350       endif
00351 !
00352 !===> ... set indices for cyclic coordinate directions
00353 !
00354          do j = 1, ndim_3d
00355          if (cyclic(j) .and. length(j) > 1) then
00356 !
00357                do n = 1, n_corners
00358 !cdir vector
00359                   do i = nbeg, nend
00360                   if (neighbors_3d (j, i, n) < grid_valid_shape(1,j)) then
00361                       neighbors_3d (j, i, n) = &
00362                       neighbors_3d (j, i, n) + length (j)
00363 
00364                   else if ( neighbors_3d (j, i, n) > grid_valid_shape(2,j)) then
00365                       neighbors_3d (j, i, n) = &
00366                       neighbors_3d (j, i, n) - length (j)
00367                   endif
00368                   end do
00369                end do ! n
00370          endif
00371          end do ! j
00372 !
00373 !===> All done
00374 !
00375 #ifdef VERBOSE
00376       print 9980, trim(ch_id)
00377 
00378       call psmile_flushstd
00379 #endif /* VERBOSE */
00380 !
00381 !  Formats:
00382 !
00383 9990  format (1x, a, ': psmile_neigh_trili_irreg2')
00384 9980  format (1x, a, ': psmile_neigh_trili_irreg2: eof')
00385 
00386 #ifdef PRISM_ASSERTION
00387 9970  format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
00388 #endif
00389 
00390       end subroutine psmile_neigh_trili_irreg2