psmile_search_nn_3d_real.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_Search_3d_real
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_search_nn_3d_real (                     &
00012                     sin_search, cos_search, z_search, distance, &
00013                     nfound, locations, n_send,                  &
00014                     x_coords, y_coords, z_coords, coords_shape, &
00015                     sin_values, cos_values, grid_valid_shape,   &
00016                     mask_array, mask_shape, mask_available,     &
00017                     tol, ierror)
00018 !
00019 ! !USES:
00020 !
00021       use PRISM
00022 !
00023       use PSMILe, dummy_interface => PSMILe_Search_nn_3d_real
00024 
00025       Implicit none
00026 
00027 !
00028 ! !DEFINED PARAMETERS:
00029 !
00030 !  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
00031 !  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
00032 !  real_earth = Earth radius in meter
00033 
00034       Integer, Parameter              :: lon = 1
00035       Integer, Parameter              :: lat = 2
00036 !
00037       Real, Parameter                 :: real_earth = 6400000.0
00038       Real, Parameter                 :: acosp1 =  1.0
00039       Real, Parameter                 :: acosm1 = -1.0
00040       Real, Parameter                 :: eps = 1d-6
00041 !
00042 ! !INPUT PARAMETERS:
00043 
00044       Integer,                     Intent (In)    :: n_send
00045 
00046 !     Number of points to be searched
00047 !
00048       Real,            Intent (In)                :: sin_search (n_send, lat)
00049 !
00050 !     Sin values of the points for which the nearest neighbours
00051 !     should be searched.
00052 !
00053       Real,            Intent (In)                :: cos_search (n_send, lat)
00054 !
00055 !     Cos values of the points for which the nearest neighbours
00056 !     should be searched.
00057 !
00058       Real,            Intent (In)                :: z_search (n_send)
00059 !
00060 !     Z values of the points for which the nearest neighbours
00061 !     should be searched.
00062 
00063       Integer,                     Intent (In)    :: coords_shape (2, ndim_3d)
00064 
00065 !     Dimension of coordinates (method) x_coords, ...
00066 
00067       Real,            Intent (In)                ::    
00068          x_coords( coords_shape(1,1):coords_shape(2,1), 
00069                    coords_shape(1,2):coords_shape(2,2), 
00070                    coords_shape(1,3):coords_shape(2,3)) 
00071       Real,            Intent (In)                ::    
00072          y_coords( coords_shape(1,1):coords_shape(2,1), 
00073                    coords_shape(1,2):coords_shape(2,2), 
00074                    coords_shape(1,3):coords_shape(2,3)) 
00075       Real,            Intent (In)                ::    
00076          z_coords( coords_shape(1,1):coords_shape(2,1), 
00077                    coords_shape(1,2):coords_shape(2,2), 
00078                    coords_shape(1,3):coords_shape(2,3)) 
00079 
00080 !     Coordinates of the method
00081 
00082       Integer,                     Intent (In)    :: mask_shape (2, ndim_3d)
00083 
00084 !     Dimension of (method) mask array
00085 
00086       Logical, Intent (In)            ::                
00087          mask_array (mask_shape (1,1):mask_shape (2,1), 
00088                      mask_shape (1,2):mask_shape (2,2), 
00089                      mask_shape (1,3):mask_shape (2,3))
00090 !     Mask of the method
00091 
00092       Logical, Intent (In)            :: mask_available
00093 
00094 !     Is mask specified in array "mask_array" ?
00095 !     mask_available = .false. : Mask is not available
00096 !                      .true.  : Mask is     available
00097 
00098       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00099 
00100 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00101 
00102       Real, Intent (In)               ::                          
00103          sin_values (grid_valid_shape(1,1):grid_valid_shape(2,1), 
00104                      grid_valid_shape(1,2):grid_valid_shape(2,2), 
00105                      grid_valid_shape(1,3):grid_valid_shape(2,3), lat)
00106 !
00107 !     Sin values of Longitudes and Latitudes (x_coords, y_coords)
00108 !
00109       Real, Intent (In)               ::                          
00110          cos_values (grid_valid_shape(1,1):grid_valid_shape(2,1), 
00111                      grid_valid_shape(1,2):grid_valid_shape(2,2), 
00112                      grid_valid_shape(1,3):grid_valid_shape(2,3), lat)
00113 !    
00114 !     Cos values of Longitudes and Latitudes (x_coords, y_coords)
00115 
00116       Real,            Intent (In)                :: tol
00117 
00118 !     Absolute tolerance for search of "identical" points
00119 !     TOL >= 0.0
00120 
00121 !
00122 ! !INPUT/OUTPUT PARAMETERS:
00123 !
00124       Real,            Intent (InOut)             :: distance (n_send)
00125 !
00126 !     Nearest neighbour distance currently found
00127 !
00128       Integer,                     Intent (InOut) :: nfound (n_send)
00129 !
00130 !     Number of nearest neighbours found per point
00131 !
00132 ! !OUTPUT PARAMETERS:
00133 !
00134       Integer,                     Intent (Out)   :: locations (ndim_3d, n_send)
00135 !
00136 !     Number of distances found per point to be searched
00137 !
00138       Integer,                     Intent (Out)   :: ierror
00139 
00140 !     Returns the error code of PSMILe_Search_nn_3d_real;
00141 !             ierror = 0 : No error
00142 !             ierror > 0 : Severe error
00143 !
00144 ! !LOCAL VARIABLES
00145 !
00146 !     ... loop variables
00147 !
00148       Integer                       :: n
00149 !
00150 !     ... for nearest neighbour control
00151 !
00152       Real                          :: real_huge, dist2, fac
00153       Integer                       :: kmin (ndim_3d)
00154 !
00155       Real, Allocatable             :: vals (:, :, :)
00156 !
00157 !     ... for error parameters
00158 !
00159       Integer, Parameter            :: nerrp = 2
00160       Integer                       :: ierrp (nerrp)
00161 !
00162 ! !DESCRIPTION:
00163 !
00164 ! Subroutine "PSMILe_Search_nn_3d_real" searches for the nearest
00165 ! point on the block boundary.
00166 !
00167 ! Subroutine "PSMILe_Search_nn_3d_real" is designed for grids of
00168 ! of type "PRISM_Irrlonlatvrt".
00169 !
00170 ! !REVISION HISTORY:
00171 !
00172 !   Date      Programmer   Description
00173 ! ----------  ----------   -----------
00174 ! 23.10.06    H. Ritzdorf  created
00175 !
00176 !EOP
00177 !----------------------------------------------------------------------
00178 !
00179 ! $Id: psmile_search_nn_3d_real.F90 2082 2009-10-21 13:31:19Z hanke $
00180 ! $Author: hanke $
00181 !
00182    Character(len=len_cvs_string), save :: mycvs = 
00183        '$Id: psmile_search_nn_3d_real.F90 2082 2009-10-21 13:31:19Z hanke $'
00184 !----------------------------------------------------------------------
00185 !
00186 !  Initialization
00187 !
00188 #ifdef VERBOSE
00189       print 9990, trim(ch_id), n_send
00190 
00191       call psmile_flushstd
00192 #endif /* VERBOSE */
00193 !
00194       ierror = 0
00195 !
00196 !===> Allocate temporary array
00197 !
00198 #if 0
00199 TODO:  allocate only data for boundary of the 3d block
00200        and test on these points
00201        whether smaller distance is available in block
00202        Currently Brute force
00203       num = 2 * ( (grid_valid_shape(2,1)-grid_valid_shape(1,1)+1)*
00204                   (grid_valid_shape(2,2)-grid_valid_shape(1,2)+1)  +
00205                   (grid_valid_shape(2,2)-grid_valid_shape(1,2)+1)*
00206                   (grid_valid_shape(2,3)-grid_valid_shape(1,3)+1)  +
00207                   (grid_valid_shape(2,3)-grid_valid_shape(1,3)+1)*
00208                   (grid_valid_shape(2,1)-grid_valid_shape(1,1)+1) )
00209 #endif
00210 
00211       Allocate (vals (grid_valid_shape(1,1):grid_valid_shape(2,1),  &
00212                       grid_valid_shape(1,2):grid_valid_shape(2,2),  &
00213                       grid_valid_shape(1,3):grid_valid_shape(2,3)), &
00214                 STAT = ierror)
00215 
00216       if (ierror > 0) then
00217          ierrp (1) = ierror
00218          ierrp (2) = (grid_valid_shape(2,1)-grid_valid_shape(1,1)+1) * &
00219                      (grid_valid_shape(2,2)-grid_valid_shape(1,2)+1)
00220 
00221          ierror = PRISM_Error_Alloc
00222          call psmile_error ( ierror, 'vals', ierrp, 2, &
00223                  __FILE__, __LINE__ )
00224          return
00225       endif
00226 !
00227 !===> Brute force
00228 !     TODO: Multigrid search
00229 !           Use psmile_neigh_nearx_... routines
00230 !
00231       real_huge = huge (distance(1))
00232 !
00233          do n = 1, n_send
00234          if (distance (n) == real_huge) then
00235             dist2 = real_huge
00236          else
00237             dist2 = distance (n) * distance (n)
00238          endif
00239 
00240          fac = real_earth + z_search (n)
00241 !
00242 !===>    ... Compute minimum distance in (lon, lat) direction
00243 !        Sehr teuer !!
00244 !
00245          vals = ( sin_values(:,:,:, lat) * sin_search(n, lat)       &
00246                 + cos_values(:,:,:, lat) * cos_search(n, lat)       &
00247                   * (cos_values(:,:,:, lon) * cos_search(n, lon)    &
00248                     +sin_values(:,:,:, lon) * sin_search(n, lon)) )
00249 !
00250 #ifdef PRISM_ASSERTION
00251          if (minval (vals) < acosm1 - eps .or. &
00252              maxval (vals) > acosp1 + eps) then
00253             print *, 'min vals', minloc (vals), minval (vals)
00254             print *, 'max vals', maxloc (vals), maxval (vals)
00255             call psmile_assert (__FILE__, __LINE__, &
00256                                 'invalid value for acos')
00257          endif
00258 #endif
00259 
00260          vals = min (acosp1, vals)
00261          vals = max (acosm1, vals)
00262 
00263          vals = fac * acos (vals)
00264          vals = vals * vals
00265 !
00266 !===> ... Add distance in Z-direction
00267 !
00268          vals = vals + &
00269             (z_coords (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00270                        grid_valid_shape(1,2):grid_valid_shape(2,2), &
00271                        grid_valid_shape(1,3):grid_valid_shape(2,3)) &
00272              - z_search (n)) ** 2
00273 !
00274          if (mask_available) then
00275             kmin = minloc (vals, &
00276                      mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00277                                  grid_valid_shape(1,2):grid_valid_shape(2,2), &
00278                                  grid_valid_shape(1,3):grid_valid_shape(2,3)) )
00279          else
00280             kmin = minloc (vals)
00281          endif
00282 !
00283          if (vals(kmin(1), kmin(2), kmin(3)) >= dist2) cycle
00284 !
00285 !===>    ... Point with smaller distance found
00286 !
00287          locations (:, n) = kmin (:)
00288 !
00289          nfound (n) = 1
00290 !
00291 #ifdef DEBUGX
00292          print *, 'psmile_search_nn_3d_real: n, old, dist2', &
00293                   n, dist2, vals(kmin(1), kmin(2), kmin(3))
00294          print *, 'psmile_search_nn_3d_real: n, loc', &
00295                   n, locations (:, n)
00296 #endif
00297 
00298          distance (n) = sqrt ( vals(kmin(1), kmin(2), kmin(3)) )
00299             
00300          enddo ! n
00301 !
00302       Deallocate (vals)
00303 !
00304 !===> All done
00305 !
00306 #ifdef VERBOSE
00307       print 9980, trim(ch_id), ierror
00308 
00309       call psmile_flushstd
00310 #endif /* VERBOSE */
00311 !
00312       return
00313 !
00314 !  Formats:
00315 !
00316 #ifdef VERBOSE
00317 9990 format (1x, a, ': psmile_search_nn_3d_real: n_send', i7)
00318 9980 format (1x, a, ': psmile_search_nn_3d_real: eof  ierror =', i4)
00319 #endif /* VERBOSE */
00320 
00321       end subroutine PSMILe_Search_nn_3d_real