psmile_neigh_trili_gauss2_irreg.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_gauss2_irreg
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_trili_gauss2_irreg (                 &
00012                grid_id, grid_valid_shape, interpolation_mode,      &
00013                srcloc, srclocz, virtual_cell, nlocs, nloc, nprev, &
00014                virtual_cell_available, neighbors_3d, num_neigh, ierror)
00015 !
00016 ! !USES:
00017 !
00018       use PRISM_constants
00019       use PSMILe, dummy_interface => PSMILe_Neigh_trili_gauss2_irreg
00020       use psmile_grid_reduced_gauss
00021 #ifdef DEBUG_TRACE
00022       use psmile_debug_trace
00023 #endif
00024 
00025       implicit none
00026 !
00027 ! !INPUT PARAMETERS:
00028 !
00029       Integer, Intent (In)            :: nloc
00030 !
00031 !     Total number of locations to be transferred
00032 !
00033       Integer, Intent (In)            :: nlocs (ndim_2d)
00034 !
00035 !     Number of locations in each direction
00036 !
00037       Integer, Intent (In)            :: srcloc  (ndim_2d, nlocs (1))
00038       Integer, Intent (In)            :: srclocz (nlocs (2))
00039 !
00040 !     Indices of the (method) grid cell in which the point was found.
00041 
00042       Integer, Intent (In)            :: nprev
00043 !
00044 !     Number of locations for which neighbours were already searched.
00045 !
00046       Integer, Intent (In)            :: virtual_cell (nlocs(1))
00047 
00048 !     Code whether destination point is located in virtual cell with
00049 !
00050 !     virtual_cell (n) = 0 : Not located in a virtual cell
00051 !     virtual_cell (n) > 0 :     Located in a virtual cell
00052 !                      = 1 : (-1,-1) direction
00053 !                      = 2 : ( 0,-1) direction
00054 !                      = 3 : ( 1,-1) direction
00055 !                      = 4 : (-1, 0) direction
00056 !                      = 6 : ( 1, 0) direction
00057 !                      = 7 : (-1, 1) direction
00058 !                      = 8 : ( 0, 1) direction
00059 !                      = 9 : ( 1, 1) direction
00060 
00061       Integer, Intent (In)            :: grid_id
00062 
00063 !     Grid handle
00064 
00065       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00066 
00067 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00068 !
00069       Integer,            Intent (In) :: interpolation_mode
00070 !
00071 !     Specifies the interpolation mode with
00072 !        interpolation_mode = 3 : Tri-linear or bilinear/linear interpolation
00073 !                                 (determines 8 corners)
00074 !        interpolation_mode = 2 : Bilinear interpolation
00075 !                                 (determines 4 corners)
00076 !        interpolation_mode = 1 : Linear   interpolation
00077 !                                 (determines 2 corners)
00078 !
00079       Integer,            Intent (In) :: num_neigh
00080 !
00081 !     Last dimension of neighbors array "neighbors_3d"
00082 !
00083       Logical,            Intent (In) :: virtual_cell_available
00084 !
00085 !     Is the virtual cell info available ?
00086 !
00087 ! !OUTPUT PARAMETERS:
00088 !
00089       Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)
00090 
00091 !     Indices of neighbor locations (interpolation bases)
00092 
00093       Integer, Intent (Out)           :: ierror
00094 
00095 !     Returns the error code of PSMILE_Neigh_trili_gauss2_irreg;
00096 !             ierror = 0 : No error
00097 !             ierror > 0 : Severe error
00098 !
00099 ! !DEFINED PARAMETERS:
00100 !
00101       integer, parameter              :: shift_3d(ndim_3d, 9) = 
00102                                              reshape ((/-1,-1,0, 0,-1,0, +1,-1,0,
00103                                                         -1, 0,0, 0, 0,0, +1, 0,0,
00104                                                         -1,+1,0, 0,+1,0, +1,+1,0/), 
00105                                                       (/ndim_3d,9/))
00106 !
00107 ! !LOCAL VARIABLES
00108 !
00109 !     ... for locations searched
00110 !
00111       Integer                         :: i, n
00112             Integer                         :: n_corners
00113 !
00114 !     ... for error handling
00115 !
00116       Integer, Parameter              :: nerrp = 2
00117       Integer                         :: ierrp (nerrp)
00118 #ifdef DEBUG_TRACE
00119 !
00120 !     ... for debugging
00121 !
00122       integer                         :: ictl_neigh_global_1d(4), ictl_neigh_local_1d(4)
00123 #endif
00124 !
00125 ! !DESCRIPTION:
00126 !
00127 ! Subroutine "psmile_neigh_trili_gauss2_irreg" generates the stencil for the cubic
00128 ! interpolation and stores the local 1d coordinates the stencil in neighbors_3d.
00129 ! In addition it stores the information whether a stencil covers the pole or
00130 ! not in neigh_bascule.
00131 !
00132 ! The subroutine is designed for (source) grids of type "PRISM_Gaussreduced_regvert"
00133 !  and target grids of type "PRISM_Irrlonlat_regvrt".
00134 !
00135 ! !REVISION HISTORY:
00136 !
00137 !   Date      Programmer   Description
00138 ! ----------  ----------   -----------
00139 ! 03.07.21    H. Ritzdorf  created
00140 ! 18.02.11    M. Hanke     rewitten based on psmile_neigh_tricu_gauss2_irreg.F90
00141 !
00142 !EOP
00143 !----------------------------------------------------------------------
00144 !
00145 !  $Id: psmile_neigh_trili_gauss2_irreg.F90 3011 2011-03-10 17:50:02Z hanke $
00146 !  $Author: hanke $
00147 !
00148    Character(len=len_cvs_string), save :: mycvs = 
00149        '$Id: psmile_neigh_trili_gauss2_irreg.F90 3011 2011-03-10 17:50:02Z hanke $'
00150 !
00151 !----------------------------------------------------------------------
00152 
00153 #ifdef VERBOSE
00154       print 9990, trim(ch_id)
00155       call psmile_flushstd
00156 #endif /* VERBOSE */
00157 
00158 !
00159 !  Initialization
00160 !
00161       ierror  = 0
00162       if (interpolation_mode == 2) then
00163          n_corners = 4
00164       else if (interpolation_mode == 3) then
00165          n_corners = 8
00166       else
00167          ierror = PRISM_Error_Internal
00168          call psmile_error ( ierror, 'unsupported interpolation mode in psmile_neigh_trili', &
00169                              (/interpolation_mode/), 1, &
00170                              __FILE__, __LINE__ )
00171       endif
00172 !
00173 !===> inital checks
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 !
00195 !cdir vector
00196          do i = 1, nlocs (1)
00197             if ( abs(srcloc(1,i)) < grid_valid_shape(1,1) .or. &
00198                  abs(srcloc(1,i)) > grid_valid_shape(2,1) .or. &
00199                  abs(srcloc(2,i)) < grid_valid_shape(1,2) .or. &
00200                  abs(srcloc(2,i)) > grid_valid_shape(2,2)) then
00201 
00202                print *, "grid valid shape is ", grid_valid_shape
00203                print *, "Incorrect index in srcloc, i", i, srcloc(1:2,i)
00204                call psmile_assert (__FILE__, __LINE__, 'wrong index')
00205             endif
00206          end do
00207 
00208          if (i <= nlocs(1)) then
00209 
00210          endif
00211 
00212 !cdir vector
00213          do i = 1, nlocs (2)
00214             if ( srclocz(i) < grid_valid_shape(1,3)-1 .or. &
00215                  srclocz(i) > grid_valid_shape(2,3)) then
00216 
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          end do
00222 #endif
00223 !
00224 !===> Initialization
00225 !
00226 !
00227       ierror = 0
00228 !
00229 !==> Actual work
00230 !
00231       ! 1st generate the 3d stencil
00232 
00233       neighbors_3d(:,:,1) = psmile_gauss_local_1d_to_3d (grid_id, abs(srcloc(1,:)), nloc)
00234 
00235       if (virtual_cell_available) then
00236          ! for all locations
00237          do i = 1, nlocs(1)
00238 
00239             n = i + nprev
00240 
00241             if (virtual_cell(i) == 0) then
00242 
00243                neighbors_3d(:,n,1:4) = psmile_gauss_get_bili_stencil (grid_id, neighbors_3d(:,n,1))
00244 
00245             else
00246 
00247                neighbors_3d(:,n,1:4) = psmile_gauss_shift_bili_stencil (grid_id, neighbors_3d(:,n,1), &
00248                                                                         shift_3d(:,virtual_cell(i)))
00249 
00250             endif
00251          enddo ! �
00252       else
00253          do i = 1, nlocs(1)
00254 
00255             n = i + nprev
00256 
00257             neighbors_3d(:,n,1:4) = psmile_gauss_get_bili_stencil (grid_id, neighbors_3d(:,n,1))
00258          enddo ! i
00259       endif
00260 
00261 #ifdef DEBUG_TRACE
00262       if (ictl > nprev .and. ictl <= nprev + nlocs(1)) then
00263          print *, "neighbors_3d (3d):"
00264          do i = 1, 4
00265             print *, ' neighbors_3d (:,ictl,i) ', i, ':', neighbors_3d (:, ictl, i)
00266          enddo
00267       endif
00268 #endif
00269 
00270       ! 2rd transform the 3d coordinates into local 1d ones
00271 
00272       do i = 1, 4
00273 
00274          neighbors_3d(1,nprev+1:nprev+nlocs(1),i) = &
00275             psmile_gauss_3d_to_global_1d (grid_id, neighbors_3d(:,nprev+1:nprev+nlocs(1),i), nlocs(1))
00276          neighbors_3d(2,nprev+1:nprev+nlocs(1),i) = 1
00277 #ifdef DEBUG_TRACE
00278          if (ictl > nprev .and. ictl <= nprev + nlocs(1)) then
00279             ictl_neigh_global_1d(i) = neighbors_3d (1, ictl, i)
00280          endif
00281 #endif
00282 
00283          neighbors_3d(1,nprev+1:nprev+nlocs(1),i) = &
00284             psmile_gauss_1d_global_to_local (grid_id, neighbors_3d(1,nprev+1:nprev+nlocs(1),i), nlocs(1), 0)
00285 
00286 #ifdef DEBUG_TRACE
00287          if (ictl > nprev .and. ictl <= nprev + nlocs(1)) then
00288             ictl_neigh_local_1d(i) = neighbors_3d (1, ictl, i)
00289          endif
00290 #endif
00291       enddo
00292 
00293 #ifdef DEBUG_TRACE
00294       if (ictl > nprev .and. ictl <= nprev + nlocs(1)) then
00295          print *, "neighbors_3d (1d global):"
00296          print *, ' neighbors_3d (1,ictl,1:4) ', ictl_neigh_global_1d
00297          print *, "neighbors_3d (1d local):"
00298          print *, ' neighbors_3d (1,ictl,1:4) ', ictl_neigh_local_1d
00299       endif
00300 #endif
00301 !
00302       if ( interpolation_mode == ndim_3d ) then
00303 !cdir vector
00304          do i = 1, nloc
00305             neighbors_3d (1:2, i, 5:8) = neighbors_3d (1:2, i, 1:4)
00306          enddo
00307       endif
00308 !
00309 !===> Duplicate first two index for third direction
00310 !
00311       do i = 2, nlocs (2)
00312          neighbors_3d (1:2, nprev+(i-1)*nlocs(1)+1:nprev+i*nlocs(1),:) = &
00313          neighbors_3d (1:2, nprev+               1:nprev+  nlocs(1),:)
00314       end do
00315 !
00316 !===> Set third index
00317 !
00318       do i = 1, nlocs (2)
00319          neighbors_3d (3, nprev+(i-1)*nlocs(1)+1:nprev+i*nlocs(1), :) = srclocz (i)
00320       end do
00321 !
00322 !===> ... K indices of neighbour indices
00323 !
00324       if (interpolation_mode == ndim_3d .and. &
00325           grid_valid_shape(2,3) - grid_valid_shape(1,3) > 0) then
00326          do n = 5, 8
00327                neighbors_3d (3, nprev+1:nprev+nlocs(1)*nlocs(2),n) = &
00328                   neighbors_3d (3, nprev+1:nprev+nlocs(1)*nlocs(2), 1) + 1
00329          enddo
00330       endif
00331 !
00332 !===> All done
00333 !
00334 #ifdef VERBOSE
00335       print 9980, trim(ch_id)
00336 
00337       call psmile_flushstd
00338 #endif /* VERBOSE */
00339 !
00340 !  Formats:
00341 !
00342 9990  format (1x, a, ': psmile_neigh_trili_gauss2_irreg')
00343 9980  format (1x, a, ': psmile_neigh_trili_gauss2_irreg: eof')
00344 
00345 #ifdef PRISM_ASSERTION
00346 9970  format (1x, a, ': number of neighbors', i2, '; expected at least', i2)
00347 #endif
00348 
00349       end subroutine psmile_neigh_trili_gauss2_irreg