psmile_search_nn_irreg2_dble.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_nn_irreg2_dble
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_search_nn_irreg2_dble (                 &
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_irreg2_dble
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 !  dble_earth = Earth radius in meter
00033 
00034       Integer, Parameter              :: lon = 1
00035       Integer, Parameter              :: lat = 2
00036 !
00037       Double Precision, Parameter     :: dble_earth = 6400000.0
00038       Double Precision, Parameter     :: acosp1 =  1.0d0
00039       Double Precision, Parameter     :: acosm1 = -1.0d0
00040       Double Precision, 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       Double Precision,            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       Double Precision,            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       Double Precision,            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       Double Precision,            Intent (In)    :: x_coords(coords_shape(1,1): 
00068                                                               coords_shape(2,1), 
00069                                                               coords_shape(1,2): 
00070                                                               coords_shape(2,2)) 
00071       Double Precision,            Intent (In)    :: y_coords(coords_shape(1,1): 
00072                                                               coords_shape(2,1), 
00073                                                               coords_shape(1,2): 
00074                                                               coords_shape(2,2)) 
00075       Double Precision,            Intent (In)    :: z_coords(coords_shape(1,3): 
00076                                                               coords_shape(2,3)) 
00077 
00078 !     Coordinates of the method
00079 
00080       Integer,                     Intent (In)    :: mask_shape (2, ndim_3d)
00081 
00082 !     Dimension of (method) mask array
00083 
00084       Logical, Intent (In)            ::                
00085          mask_array (mask_shape (1,1):mask_shape (2,1), 
00086                      mask_shape (1,2):mask_shape (2,2), 
00087                      mask_shape (1,3):mask_shape (2,3))
00088 !     Mask of the method
00089 
00090       Logical, Intent (In)            :: mask_available
00091 
00092 !     Is mask specified in array "mask_array" ?
00093 !     mask_available = .false. : Mask is not available
00094 !                      .true.  : Mask is     available
00095 
00096       Integer, Intent (In)            :: grid_valid_shape (2, ndim_3d)
00097 
00098 !     Specifies the valid block shape for the "ndim_3d"-dimensional block
00099 
00100       Double Precision, Intent (In)   :: sin_values (grid_valid_shape(1,1): 
00101                                                      grid_valid_shape(2,1), 
00102                                                      grid_valid_shape(1,2): 
00103                                                      grid_valid_shape(2,2), 
00104                                                      lat)
00105 !
00106 !     Sin values of Longitudes and Latitudes (x_coords, y_coords)
00107 !
00108       Double Precision, Intent (In)   :: cos_values (grid_valid_shape(1,1): 
00109                                                      grid_valid_shape(2,1), 
00110                                                      grid_valid_shape(1,2): 
00111                                                      grid_valid_shape(2,2), 
00112                                                      lat)
00113 !    
00114 !     Cos values of Longitudes and Latitudes (x_coords, y_coords)
00115 
00116       Double Precision,            Intent (In)    :: tol
00117 
00118 !     Absolute tolerance for search of "identical" points
00119 !     TOL >= 0.0
00120 
00121 !
00122 ! !INPUT/OUTPUT PARAMETERS:
00123 !
00124       Double Precision,            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_irreg2_dble;
00141 !             ierror = 0 : No error
00142 !             ierror > 0 : Severe error
00143 !
00144 ! !LOCAL VARIABLES
00145 !
00146 !     ... loop variables
00147 !
00148       Integer                       :: n, k
00149 !
00150 !     ... for nearest neighbour control
00151 !
00152       Double Precision              :: dble_huge, distk, dist2, fac
00153       Integer                       :: ijmin (lat), kmin (1), k1
00154 !
00155       Double Precision, Allocatable :: vals (:, :)
00156 !
00157 !     ... for masks
00158 !
00159       Double Precision              :: dist_mask, dist2_mask
00160       Integer                       :: kex
00161 !
00162       Integer                       :: ij (lat)
00163       Logical,          Allocatable :: kmask (:)
00164 !
00165 !     ... for error parameters
00166 !
00167       Integer, Parameter            :: nerrp = 2
00168       Integer                       :: ierrp (nerrp)
00169 !
00170 ! !DESCRIPTION:
00171 !
00172 ! Subroutine "PSMILe_Search_nn_irreg2_dble" searches for the nearest
00173 ! point on the block boundary.
00174 !
00175 ! Subroutine "PSMILe_Search_nn_irreg2_dble" is designed for grids of
00176 ! of type "PRISM_Irrlonlat_regvrt".
00177 !
00178 ! !REVISION HISTORY:
00179 !
00180 !   Date      Programmer   Description
00181 ! ----------  ----------   -----------
00182 ! 23.10.06    H. Ritzdorf  created
00183 !
00184 !EOP
00185 !----------------------------------------------------------------------
00186 !
00187 ! $Id: psmile_search_nn_irreg2_dble.F90 2082 2009-10-21 13:31:19Z hanke $
00188 ! $Author: hanke $
00189 !
00190    Character(len=len_cvs_string), save :: mycvs = 
00191        '$Id: psmile_search_nn_irreg2_dble.F90 2082 2009-10-21 13:31:19Z hanke $'
00192 !----------------------------------------------------------------------
00193 !
00194 !  Initialization
00195 !
00196 #ifdef VERBOSE
00197       print 9990, trim(ch_id), n_send
00198 
00199       call psmile_flushstd
00200 #endif /* VERBOSE */
00201 !
00202       ierror = 0
00203 !
00204 !===> Allocate temporary array
00205 !
00206       Allocate (vals (grid_valid_shape(1,1):grid_valid_shape(2,1),  &
00207                       grid_valid_shape(1,2):grid_valid_shape(2,2)), &
00208                 STAT = ierror)
00209 
00210       if (ierror > 0) then
00211          ierrp (1) = ierror
00212          ierrp (2) = (grid_valid_shape(2,1)-grid_valid_shape(1,1)+1) * &
00213                      (grid_valid_shape(2,2)-grid_valid_shape(1,2)+1)
00214 
00215          ierror = PRISM_Error_Alloc
00216          call psmile_error ( ierror, 'vals', ierrp, 2, &
00217                  __FILE__, __LINE__ )
00218          return
00219       endif
00220 !
00221       if (mask_available) then
00222          Allocate (kmask (grid_valid_shape(1,3):grid_valid_shape(2,3)), &
00223                    STAT = ierror)
00224 
00225          if (ierror > 0) then
00226             ierrp (1) = ierror
00227             ierrp (2) = grid_valid_shape(2,3) - grid_valid_shape(1,3) + 1
00228 
00229             ierror = PRISM_Error_Alloc
00230             call psmile_error ( ierror, 'kmask', ierrp, 2, &
00231                  __FILE__, __LINE__ )
00232             return
00233          endif
00234 
00235             do k = grid_valid_shape(1,3), grid_valid_shape(2,3)
00236             kmask (k) = any (                                              &
00237                mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1),    &
00238                            grid_valid_shape(1,2):grid_valid_shape(2,2), k) )
00239             end do
00240 
00241 #ifdef PRISM_ASSERTION
00242             if (.not. ANY(kmask)) then
00243                print *, 'kmask ', kmask
00244                call psmile_assert (__FILE__, __LINE__, &
00245                        'at least one entry of kmask should be true')
00246             endif
00247 #endif
00248       endif
00249 !
00250 !===> Brute force
00251 !     TODO: Multigrid search
00252 !
00253       dble_huge = huge (distance(1))
00254       k1 = grid_valid_shape (1,3) - 1
00255 !
00256          do n = 1, n_send
00257          if (distance (n) == dble_huge) then
00258             dist2 = dble_huge
00259          else
00260             dist2 = distance (n) * distance (n)
00261          endif
00262 
00263          fac = dble_earth + z_search (n)
00264 !
00265 !===>    ... Compute minimum distance in z-direction
00266 !
00267          if (mask_available) then
00268             kmin = minloc (abs(z_coords (grid_valid_shape(1,3):                  &
00269                                          grid_valid_shape(2,3))) - z_search (n), &
00270                            kmask)
00271          else
00272             kmin = minloc (abs(z_coords (grid_valid_shape(1,3):                  &
00273                                          grid_valid_shape(2,3))) - z_search (n))
00274          endif
00275 !
00276          kmin (1) = kmin (1) + k1
00277 !
00278          distk = abs (z_coords (kmin(1)) - z_search (n))
00279          if (distk >= distance (n)) cycle
00280 !
00281 !===>    ... Compute minimum distance in (lon, lat) direction
00282 !
00283          vals = ( sin_values(:,:, lat) * sin_search(n, lat)       &
00284                 + cos_values(:,:, lat) * cos_search(n, lat)       &
00285                   * (cos_values(:,:, lon) * cos_search(n, lon)    &
00286                     +sin_values(:,:, lon) * sin_search(n, lon)) )
00287 !
00288 #ifdef PRISM_ASSERTION
00289          if (minval (vals) < acosm1 - eps .or. &
00290              maxval (vals) > acosp1 + eps) then
00291             print *, 'min vals', minloc (vals), minval (vals)
00292             print *, 'max vals', maxloc (vals), maxval (vals)
00293             call psmile_assert (__FILE__, __LINE__, &
00294                                 'invalid value for acos')
00295          endif
00296 #endif
00297 
00298          vals = min (acosp1, vals)
00299          vals = max (acosm1, vals)
00300 
00301          vals = fac * acos (vals)
00302          vals = vals * vals + distk * distk
00303 !
00304          if (mask_available) then
00305             dist2_mask = dist2
00306             dist_mask  = distance (n)
00307 !
00308 !===> Look for an valid entry
00309 !
00310             if (kmask (kmin(1))) then
00311                ijmin = minloc (vals,                                            &
00312                        mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00313                                    grid_valid_shape(1,2):grid_valid_shape(2,2), &
00314                                    kmin (1)) )
00315                ijmin = ijmin + grid_valid_shape (1, 1:lat) - 1
00316                if (vals (ijmin(1), ijmin(2)) < dist2_mask) then
00317                   dist2_mask = vals (ijmin(1), ijmin(2))
00318                   dist_mask  = sqrt (dist_mask)
00319                endif
00320 
00321                kex = kmin (1) ! exclude kmin(1); already computed
00322                kmask (kex) = .false.
00323             else
00324                kex = grid_valid_shape(1,3) - 1
00325             endif
00326 !
00327 !
00328                do k = grid_valid_shape(1,3), grid_valid_shape(2,3)
00329                if (.not. kmask(k)) cycle
00330                if (abs (z_coords (k) - z_search (n)) >= dist_mask) cycle
00331 
00332                   ij = minloc (vals,                                          &
00333                      mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00334                                  grid_valid_shape(1,2):grid_valid_shape(2,2), &
00335                                  k) )
00336                   ij = ij + grid_valid_shape (1, 1:lat) - 1
00337 
00338                   if (vals (ij(1),ij(2)) < dist_mask) then
00339                      dist2_mask = vals (ijmin(1),ijmin(2))
00340                      dist_mask  = sqrt (dist2_mask)
00341                      ijmin = ij
00342                      kmin (1) = k
00343                   endif
00344                end do
00345 
00346             if (kex /= grid_valid_shape(1,3)-1) kmask (kex) = .true.
00347 
00348          else
00349 
00350             ijmin = minloc (vals)
00351             ijmin = ijmin + grid_valid_shape (1, 1:lat) - 1
00352 
00353          endif
00354 
00355 !
00356          if (vals(ijmin(1), ijmin(2)) >= dist2) cycle
00357 !
00358 !===>    ... Point with smaller distance found
00359 !
00360          locations (1, n) = ijmin (1)
00361          locations (2, n) = ijmin (2)
00362          locations (3, n) = kmin (1)
00363 !
00364          nfound (n) = 1
00365 !
00366 #ifdef DEBUGX
00367          print *, 'psmile_search_nn_irreg2_dble: n, old, dist2', &
00368                   n, dist2, vals (ijmin(1), ijmin(2))
00369          print *, 'psmile_search_nn_irreg2_dble: n, loc', &
00370                   n, locations (:, n)
00371 #endif
00372 
00373          distance (n) = sqrt (vals (ijmin(1), ijmin(2)))
00374             
00375          enddo ! n
00376 !
00377       Deallocate (vals)
00378       if (mask_available) Deallocate (kmask)
00379 !
00380 !===> All done
00381 !
00382 #ifdef VERBOSE
00383       print 9980, trim(ch_id), ierror
00384 
00385       call psmile_flushstd
00386 #endif /* VERBOSE */
00387 !
00388       return
00389 !
00390 !  Formats:
00391 !
00392 #ifdef VERBOSE
00393 9990 format (1x, a, ': psmile_search_nn_irreg2_dble: n_send', i7)
00394 9980 format (1x, a, ': psmile_search_nn_irreg2_dble: eof  ierror =', i4)
00395 #endif /* VERBOSE */
00396 
00397       end subroutine PSMILe_Search_nn_irreg2_dble