psmile_neigh_nearx_3d_irr3_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2008-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_nearx_3d_irr3_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_neigh_nearx_3d_irr3_dble (         &
00012                            grid_id,                        &
00013                            coords1, coords2, coords3,      &
00014                            x_coords, y_coords, z_coords,   &
00015                            coords_shape,                   &
00016                            mask_array, mask_shape,  mask_available,    &
00017                            sin_values_lon, cos_values_lon, &
00018                            sin_values_lat, cos_values_lat, &
00019                            grid_valid_shape, search_mode,  &
00020                            srcloc, nsearch,                &
00021                            nloc, nprev,                    &
00022                            neighbors_3d, num_neigh,        &
00023                            extra_search, ierror)
00024 !
00025 ! !USES:
00026 !
00027       use PRISM_constants
00028 !
00029       use PSMILe, dummy_interface => PSMILe_Neigh_nearx_3d_irr3_dble
00030 
00031       implicit none
00032 !
00033 ! !INPUT PARAMETERS:
00034 !
00035       Integer, Intent (In)            :: grid_id
00036 
00037 !     Info on the component in which the donor cells
00038 !     should be searched.
00039 
00040       Integer, Intent (In)            :: nloc
00041 !
00042 !     Total number of locations to be transferred
00043 !
00044       Integer, Intent (In)            :: nprev
00045 !
00046 !     Number of locations for which neighbours were already searched.
00047 !
00048       Integer, Intent (In)            :: nsearch
00049 !
00050 !     Number of locations for which neighbours should be searched in this
00051 !     sweep.
00052 !
00053       Integer, Intent (In)            :: srcloc (ndim_3d, nloc)
00054 
00055 !     Indices of the grid cell in which the point was found.
00056 !
00057       Double Precision, Intent (In)   :: coords1 (nloc)
00058       Double Precision, Intent (In)   :: coords2 (nloc)
00059       Double Precision, Intent (In)   :: coords3 (nloc)
00060 
00061 !     Coordinates of the target points (which were searched and found)
00062 
00063       Integer, Intent (In)            :: coords_shape (2, ndim_3d)
00064 
00065 !     Dimension of coordinates (method) x_coords, ...
00066 
00067       Double Precision, Intent (In)   :: x_coords(coords_shape(1,1): 
00068                                                   coords_shape(2,1))
00069       Double Precision, Intent (In)   :: y_coords(coords_shape(1,2): 
00070                                                   coords_shape(2,2))
00071       Double Precision, Intent (In)   :: z_coords(coords_shape(1,3): 
00072                                                   coords_shape(2,3))
00073 
00074 !     Coordinates of the method
00075 
00076       Integer, Intent (In)            :: mask_shape (2, ndim_3d)
00077 
00078 !     Dimension of (method) mask array
00079 
00080       Logical, Intent (In)            :: mask_array (mask_shape (1,1): 
00081                                                      mask_shape (2,1), 
00082                                                      mask_shape (1,2): 
00083                                                      mask_shape (2,2), 
00084                                                      mask_shape (1,3): 
00085                                                      mask_shape (2,3))
00086 !     Mask of the method
00087 
00088       Logical, Intent (In)            :: mask_available
00089 
00090 !     Is mask specified in array "mask_array" ?
00091 !     mask_available = .false. : Mask is not available
00092 !                      .true.  : Mask is     available
00093 
00094       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00095 
00096 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00097 
00098       Double Precision, Intent (In)   :: sin_values_lon (grid_valid_shape(1,1): 
00099                                                          grid_valid_shape(2,1))
00100       Double Precision, Intent (In)   :: sin_values_lat (grid_valid_shape(1,2): 
00101                                                          grid_valid_shape(2,2))
00102 !
00103 !     Sin values of Longitudes and Latitudes (x_coords, y_coords)
00104 !
00105       Double Precision, Intent (In)   :: cos_values_lon (grid_valid_shape(1,1): 
00106                                                          grid_valid_shape(2,1))
00107       Double Precision, Intent (In)   :: cos_values_lat (grid_valid_shape(1,2): 
00108                                                          grid_valid_shape(2,2))
00109 !    
00110 !     Cos values of Longitudes and Latitudes (x_coords, y_coords)
00111 !
00112       Integer, Intent (In)            :: search_mode
00113 !
00114 !     Specifies the search mode for nearest neigbours with
00115 !        search_mode = 3 : Full 3d-search
00116 !        search_mode = 2 : Search in 2d-hyperplane (1st 2 direction lon, lat)
00117 !        search_mode = 1 : Search in 1d-plane
00118 !
00119       Integer, Intent (In)            :: num_neigh
00120 !
00121 !     Last dimension of neighbors array "neighbors_3d" and
00122 !     number of neighbors to be searched.
00123 !
00124 ! !INPUT/OUTPUT PARAMETERS:
00125 !
00126       Type (Extra_search_info)        :: extra_search
00127 !
00128 !     Number of locations where
00129 !       (*) required mask values were not true
00130 !
00131 ! !OUTPUT PARAMETERS:
00132 !
00133       Integer, Intent (Out)           :: neighbors_3d (ndim_3d, nloc, num_neigh)
00134 
00135 !     Indices of neighbor locations (interpolation bases)
00136 !
00137       Integer, Intent (Out)           :: ierror
00138 
00139 !     Returns the error code of PSMILE_Neigh_nearx_3d_irr3_dble;
00140 !             ierror = 0 : No error
00141 !             ierror > 0 : Severe error
00142 !
00143 ! !DEFINED PARAMETERS:
00144 !
00145 !  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
00146 !  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
00147 !
00148       Integer, Parameter              :: lon = 1
00149       Integer, Parameter              :: lat = 2
00150 !
00151 ! !LOCAL VARIABLES
00152 !
00153 !     ... For locations controlled
00154 !
00155       Double Precision, Pointer       :: sin_search (:, :)
00156       Double Precision, Pointer       :: cos_search (:, :)
00157       Double Precision, Pointer       :: z_search (:)
00158 !
00159 !     ... (1d) Indices in (entire) data to be extra searched
00160 !
00161 #ifdef PRISM_ASSERTION
00162       Integer                         :: i
00163 #endif
00164       Integer                         :: j
00165       Integer                         :: ibeg, iend
00166       Integer                         :: jbeg, jend
00167 !
00168       Integer, Pointer                :: indices (:)
00169 !
00170 !     ... 3d Indices in local data to be extra searched
00171 !
00172       Integer                         :: ijk0 (ndim_3d)
00173       Integer, Pointer                :: ijk (:, :)
00174       Integer                         :: width (ndim_3d)
00175 !
00176 !     ... for error handling
00177 !
00178       Integer, Parameter              :: nerrp = 2
00179       Integer                         :: ierrp (nerrp)
00180 !
00181 !
00182 ! !DESCRIPTION:
00183 !
00184 ! Subroutine "PSMILe_Neigh_nearx_3d_irr3_dble" searches the next "num_neigh"
00185 ! nearest neighbours on the method-grid (x_coords, y_coords, z_coords)
00186 ! for the extra indices specified in "extra_search%indices".
00187 !
00188 ! There are 3 different of search modes supported by the routine:
00189 ! search_mode = 3: The nearest neighbours are searched within the
00190 !                  the 3d region.
00191 ! search_mode = 2: The nearest neighbours are searched only in the
00192 !                  2d hyperplane (lon- and lat-direction).
00193 ! search_mode = 1: The nearest neighbours are searched only in the
00194 !                  1d plane.
00195 !
00196 ! Subroutine "PSMILe_Neigh_nearx_3d_irr3_dble" is designed
00197 ! (*) for source grids of type "PRISM_Reglonlatvrt" and
00198 ! (*) source locations "srcloc" which are stored as a full 3d vector
00199 !     containing the 3d indices of the cells found for the target points
00200 !     "coords" to be searched.
00201 !
00202 ! !REVISION HISTORY:
00203 !
00204 !   Date      Programmer   Description
00205 ! ----------  ----------   -----------
00206 ! 03.07.21    H. Ritzdorf  created
00207 !
00208 !EOP
00209 !----------------------------------------------------------------------
00210 !
00211 !  $Id: psmile_neigh_nearx_3d_irr3_dble.F90 2966 2011-02-18 09:47:30Z hanke $
00212 !  $Author: hanke $
00213 !
00214    Character(len=len_cvs_string), save :: mycvs = 
00215        '$Id: psmile_neigh_nearx_3d_irr3_dble.F90 2966 2011-02-18 09:47:30Z hanke $'
00216 !                            
00217 !----------------------------------------------------------------------
00218 !
00219 !  Initialization
00220 !
00221 #ifdef VERBOSE
00222       print 9990, trim(ch_id)
00223 
00224       call psmile_flushstd
00225 #endif /* VERBOSE */
00226 !
00227       ierror = 0
00228 !
00229       ibeg = nprev + 1
00230       iend = nprev + nsearch
00231 !
00232       indices => extra_search%indices
00233 !
00234 !===> Get range in indices to be controlled
00235 !     Note: indices are locations in the full 3d cube and
00236 !           contains the previous numbers
00237 !           TODO: DAS SOLLTE in lokale Indices aendern
00238 !
00239          do jbeg = 1, extra_search%n_extra
00240          if (indices(jbeg) >= ibeg) exit
00241          end do
00242 !
00243       if (jbeg > extra_search%n_extra) then
00244 #ifdef VERBOSE
00245          print 9980, trim(ch_id), ierror
00246          call psmile_flushstd
00247 #endif /* VERBOSE */
00248          return
00249       endif
00250 !
00251          do jend = jbeg, extra_search%n_extra
00252          if (indices(jend) > iend) exit
00253          end do
00254       jend = jend - 1
00255 !
00256       if (jbeg > jend) then
00257 #ifdef VERBOSE
00258          print 9980, trim(ch_id), ierror
00259          call psmile_flushstd
00260 #endif /* VERBOSE */
00261          return
00262       endif
00263 !
00264       ibeg = indices(jbeg)
00265       iend = indices(jend)
00266 !
00267 #ifdef PRISM_ASSERTION
00268 !
00269 !===> Internal control
00270 !
00271       if (search_mode < 1 .or. search_mode > ndim_3d) then
00272          print *, trim(ch_id), " search_mode ", search_mode
00273          call psmile_assert (__FILE__, __LINE__, &
00274                              'search_mode is out of range ')
00275       endif
00276 !
00277       if (nloc < nprev+nsearch ) then
00278          print *, trim(ch_id), " nloc, nprev, nsearch ", nloc, nprev, nsearch
00279          call psmile_assert (__FILE__, __LINE__, &
00280                              'nloc < nprev + nsearch ')
00281       endif
00282 !
00283 !===> Are the locations within the correct shape ?
00284 !
00285 !cdir vector
00286          do j = jbeg, jend
00287          i = indices(j)
00288 !
00289          if (srcloc(1,i) < coords_shape(1,1) - 1 .or. &
00290              srcloc(1,i) > coords_shape(2,1)     .or. &
00291              srcloc(2,i) < coords_shape(1,2) - 1 .or. &
00292              srcloc(2,i) > coords_shape(2,2)     .or. &
00293              srcloc(3,i) < coords_shape(1,3) - 1 .or. &
00294              srcloc(3,i) > coords_shape(2,3)) exit
00295             end do
00296 !
00297       if (j <= jend) then
00298          print *, "i, j, jbeg, jend", i, j, jbeg, jend
00299          print *, "srcloc ", srcloc(:, i)
00300          print *, "shape  ", coords_shape
00301             call psmile_assert (__FILE__, __LINE__, &
00302                                 'wrong index in srcloc')
00303          endif
00304 #endif
00305 !
00306 !-----------------------------------------------------------------------
00307 !     Allocate and set temporary array which is used to transform
00308 !     the coordinates of points to be searched.
00309 !     (*) Longitudes and Latitudes are transformed
00310 !         from degrees into radients.
00311 !     (*) The z-values are currently not transformed
00312 !         ??? Whats about PRISM_Irrlonlat_sigmavrt, ...
00313 !
00314 !     Is it possible to save the values, in order to be reused.
00315 !-----------------------------------------------------------------------
00316 !
00317       Allocate (sin_search(jbeg:jend, lat), STAT = ierror)
00318 
00319       if ( ierror > 0 ) then
00320          ierrp (1) = ierror
00321          ierrp (2) = (jend-jbeg+1) * lat
00322 
00323          ierror = PRISM_Error_Alloc
00324          call psmile_error ( ierror, 'sin_search', &
00325                              ierrp, 2, __FILE__, __LINE__ )
00326          return
00327       endif
00328 
00329       Allocate (cos_search(jbeg:jend, lat), STAT = ierror)
00330 
00331       if ( ierror > 0 ) then
00332          ierrp (1) = ierror
00333          ierrp (2) = (jend-jbeg+1) * lat
00334 
00335          ierror = PRISM_Error_Alloc
00336          call psmile_error ( ierror, 'cos_search', &
00337                              ierrp, 2, __FILE__, __LINE__ )
00338          return
00339       endif
00340 !
00341       Allocate (z_search(jbeg:jend), STAT = ierror)
00342 
00343       if ( ierror > 0 ) then
00344          ierrp (1) = ierror
00345          ierrp (2) = (jend-jbeg+1)
00346 
00347          ierror = PRISM_Error_Alloc
00348          call psmile_error ( ierror, 'z_search', &
00349                              ierrp, 2, __FILE__, __LINE__ )
00350          return
00351       endif
00352 !
00353       Allocate (ijk(ndim_3d, jbeg:jend), STAT = ierror)
00354 
00355       if ( ierror > 0 ) then
00356          ierrp (1) = ierror
00357          ierrp (2) = (jend-jbeg+1) * ndim_3d
00358 
00359          ierror = PRISM_Error_Alloc
00360          call psmile_error ( ierror, 'ijk', &
00361                              ierrp, 2, __FILE__, __LINE__ )
00362          return
00363       endif
00364 !
00365 !===> Compute values required to compute distances
00366 !     Note: In the previous version (i.e. without using
00367 !           subroutine PSMILE_Neigh_nearx_sub_reg_dble),
00368 !           z_search wasn't allocated and coords3 was used
00369 !           within the body in order to compute the distance.
00370 !           PSMILE_Neigh_nearx_sub_reg_dble was included to 
00371 !           simplify the code.
00372 !
00373 !cdir vector
00374          do j = jbeg, jend
00375          sin_search (j, lon) = coords1 (indices(j)-nprev) * dble_deg2rad
00376          end do
00377 !
00378 !cdir vector
00379          do j = jbeg, jend
00380          sin_search (j, lat) = coords2 (indices(j)-nprev) * dble_deg2rad
00381          end do
00382 !
00383       cos_search = cos (sin_search)
00384       sin_search = sin (sin_search)
00385 !
00386 !cdir vector
00387          do j = jbeg, jend
00388          z_search (j) = coords3 (indices(j)-nprev)
00389          end do
00390 !
00391 !-----------------------------------------------------------------------
00392 !     Compute coarse grid indices
00393 !-----------------------------------------------------------------------
00394 !
00395       ijk0 = Grids(grid_id)%ijk0
00396 !
00397 !===> Get coarse grid index on level nlev-2.
00398 !     By coarsening with a factor of 4
00399 !
00400       width (1:search_mode) = 3
00401       width (search_mode+1:ndim_3d) = 0
00402 !
00403       ijk (1, jbeg:jend) = ijk0(1) + ((srcloc(1, indices(jbeg:jend)) - &
00404                                        ijk0(1)) / 4) * 4
00405 !
00406       if (search_mode >= 2) then
00407          ijk (2, jbeg:jend) = ijk0(2) + ((srcloc(2, indices(jbeg:jend)) - &
00408                                          ijk0(2)) / 4) * 4
00409       else
00410          ijk (2, jbeg:jend) = srcloc(2, indices(jbeg:jend))
00411       endif
00412 !
00413       if (search_mode >= 3) then
00414          ijk (3, jbeg:jend) = ijk0(3) + ((srcloc(3, indices(jbeg:jend)) - &
00415                                          ijk0(3)) / 4) * 4
00416       else
00417          ijk (3, jbeg:jend) = srcloc(3, indices(jbeg:jend))
00418       endif
00419 !
00420 !===> Compute nearest neigbours for all points
00421 !     specified in "indices(jbeg:jend)"
00422 !
00423       call psmile_neigh_nearx_sub_reg_dble (               &
00424                            grid_id,                        &
00425                            x_coords, y_coords, z_coords,   &
00426                            coords_shape,                   &
00427                            mask_array, mask_shape,  mask_available,    &
00428                            sin_values_lon, cos_values_lon, &
00429                            sin_values_lat, cos_values_lat, &
00430                            grid_valid_shape, search_mode,  &
00431                            neighbors_3d, num_neigh, nloc,  &
00432                            extra_search,                   &
00433                            ijk, sin_search, cos_search,    &
00434                            z_search, jbeg, jend, width, ierror)
00435 !
00436 !===> Deallocate memory
00437 !
00438       Deallocate (ijk)
00439       Deallocate (z_search)
00440       Deallocate (cos_search)
00441       Deallocate (sin_search)
00442 !
00443 !===> All done
00444 !
00445 #ifdef VERBOSE
00446       print 9980, trim(ch_id), ierror
00447 
00448       call psmile_flushstd
00449 #endif /* VERBOSE */
00450 !
00451 !  Formats:
00452 !
00453 9990 format (1x, a, ': psmile_neigh_nearx_3d_irr3_dble')
00454 9980 format (1x, a, ': psmile_neigh_nearx_3d_irr3_dble: eof, ierror =', i3)
00455 
00456       end subroutine PSMILe_Neigh_nearx_3d_irr3_dble