psmile_search_nn_3d_reg_real.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2007-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_3d_reg_real
00008 !
00009 ! !INTERFACE:
00010 
00011       subroutine psmile_search_nn_3d_reg_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_lon, cos_values_lon,             &
00016                     sin_values_lat, cos_values_lat,             &
00017                     grid_valid_shape,                           &
00018                     mask_array, mask_shape, mask_available,     &
00019                     tol, ierror)
00020 !
00021 ! !USES:
00022 !
00023       use PRISM
00024 !
00025       use PSMILe, dummy_interface => PSMILe_Search_nn_3d_reg_real
00026 
00027       Implicit none
00028 
00029 !
00030 ! !DEFINED PARAMETERS:
00031 !
00032 !  lon   = Index of Longitudes in arrays "sin_values" and "cos_values"
00033 !  lat   = Index of Latitudes  in arrays "sin_values" and "cos_values"
00034 !  real_earth = Earth radius in meter
00035 
00036       Integer, Parameter              :: lon = 1
00037       Integer, Parameter              :: lat = 2
00038 !
00039       Real, Parameter                 :: real_earth = 6400000.0
00040       Real, Parameter                 :: acosp1 =  1.0
00041       Real, Parameter                 :: acosm1 = -1.0
00042       Real, Parameter                 :: eps = 1d-6
00043 !
00044 ! !INPUT PARAMETERS:
00045 
00046       Integer,                     Intent (In)    :: n_send
00047 
00048 !     Number of points to be searched
00049 !
00050       Real,            Intent (In)                :: sin_search (n_send, lat)
00051 !
00052 !     Sin values of the points for which the nearest neighbours
00053 !     should be searched.
00054 !
00055       Real,            Intent (In)                :: cos_search (n_send, lat)
00056 !
00057 !     Cos values of the points for which the nearest neighbours
00058 !     should be searched.
00059 !
00060       Real,            Intent (In)                :: z_search (n_send)
00061 !
00062 !     Z values of the points for which the nearest neighbours
00063 !     should be searched.
00064 
00065       Integer,                     Intent (In)    :: coords_shape (2, ndim_3d)
00066 
00067 !     Dimension of coordinates (method) x_coords, ...
00068 
00069       Real,            Intent (In)                ::    
00070          x_coords (coords_shape(1,1):coords_shape(2,1))
00071       Real,            Intent (In)                ::    
00072          y_coords (coords_shape(1,2):coords_shape(2,2))
00073       Real,            Intent (In)                ::    
00074          z_coords (coords_shape(1,3):coords_shape(2,3)) 
00075 
00076 !     Coordinates of the method
00077 
00078       Integer,                     Intent (In)    :: mask_shape (2, ndim_3d)
00079 
00080 !     Dimension of (method) mask array
00081 
00082       Logical, Intent (In)            ::                
00083          mask_array (mask_shape (1,1):mask_shape (2,1), 
00084                      mask_shape (1,2):mask_shape (2,2), 
00085                      mask_shape (1,3):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 
00099       Real, Intent (In)               :: 
00100          sin_values_lon (grid_valid_shape(1,1):grid_valid_shape(2,1))
00101       Real, Intent (In)               :: 
00102          sin_values_lat (grid_valid_shape(1,2):grid_valid_shape(2,2))
00103 !
00104 !     Sin values of Longitudes and Latitudes (x_coords, y_coords)
00105 !
00106       Real, Intent (In)               :: 
00107          cos_values_lon (grid_valid_shape(1,1):grid_valid_shape(2,1))
00108       Real, Intent (In)               :: 
00109          cos_values_lat (grid_valid_shape(1,2):grid_valid_shape(2,2))
00110 !    
00111 !     Cos values of Longitudes and Latitudes (x_coords, y_coords)
00112 
00113       Real,            Intent (In)                :: tol
00114 
00115 !     Absolute tolerance for search of "identical" points
00116 !     TOL >= 0.0
00117 
00118 !
00119 ! !INPUT/OUTPUT PARAMETERS:
00120 !
00121       Real,            Intent (InOut)             :: distance (n_send)
00122 !
00123 !     Nearest neighbour distance currently found
00124 !
00125       Integer,                     Intent (InOut) :: nfound (n_send)
00126 !
00127 !     Number of nearest neighbours found per point
00128 !
00129 ! !OUTPUT PARAMETERS:
00130 !
00131       Integer,                     Intent (Out)   :: locations (ndim_3d, n_send)
00132 !
00133 !     Number of distances found per point to be searched
00134 !
00135       Integer,                     Intent (Out)   :: ierror
00136 
00137 !     Returns the error code of PSMILe_Search_nn_3d_reg_real;
00138 !             ierror = 0 : No error
00139 !             ierror > 0 : Severe error
00140 !
00141 ! !LOCAL VARIABLES
00142 !
00143 !     ... loop variables
00144 !
00145       Integer                       :: j, k, n
00146       Integer                       :: jend
00147 !
00148 !     ... for nearest neighbour control
00149 !
00150       Real                          :: real_huge, distk, dist2, fac
00151       Integer                       :: ijmin (lat), kmin (1), k1
00152 !
00153       Real, Allocatable             :: vals (:, :)
00154 !
00155 !     ... for masks
00156 !
00157       Real                          :: dist_mask, dist2_mask
00158       Integer                       :: kex
00159 !
00160       Integer                       :: ij (lat)
00161       Logical,          Allocatable :: kmask (:)
00162 !
00163 !     ... for error parameters
00164 !
00165       Integer, Parameter            :: nerrp = 2
00166       Integer                       :: ierrp (nerrp)
00167 !
00168 ! !DESCRIPTION:
00169 !
00170 ! Subroutine "PSMILe_Search_nn_3d_reg_real" searches for the nearest
00171 ! point on the block boundary.
00172 !
00173 ! Subroutine "PSMILe_Search_nn_3d_reg_real" is designed for grids of
00174 ! of type "PRISM_Reglonlatvrt".
00175 !
00176 ! !REVISION HISTORY:
00177 !
00178 !   Date      Programmer   Description
00179 ! ----------  ----------   -----------
00180 ! 23.10.06    H. Ritzdorf  created
00181 !
00182 !EOP
00183 !----------------------------------------------------------------------
00184 !
00185 ! $Id: psmile_search_nn_3d_reg_real.F90 2082 2009-10-21 13:31:19Z hanke $
00186 ! $Author: hanke $
00187 !
00188    Character(len=len_cvs_string), save :: mycvs = 
00189        '$Id: psmile_search_nn_3d_reg_real.F90 2082 2009-10-21 13:31:19Z hanke $'
00190 !----------------------------------------------------------------------
00191 !
00192 !  Initialization
00193 !
00194 #ifdef VERBOSE
00195       print 9990, trim(ch_id), n_send
00196 
00197       call psmile_flushstd
00198 #endif /* VERBOSE */
00199 ! TODO --- kmask alle false ?
00200 !
00201       ierror = 0
00202 !
00203 #ifdef PRISM_ASSERTION
00204       if (grid_valid_shape(1,2) > grid_valid_shape(2,2)) then
00205          print *, 'grid_valid_shape', grid_valid_shape
00206          call psmile_assert (__FILE__, __LINE__, &
00207                  "grid shape shouldn't be empty")
00208       endif
00209 #endif
00210 !
00211 !===> Allocate temporary array
00212 !
00213       Allocate (vals (grid_valid_shape(1,1):grid_valid_shape(2,1),  &
00214                       grid_valid_shape(1,2):grid_valid_shape(2,2)), &
00215                 STAT = ierror)
00216 
00217       if (ierror > 0) then
00218          ierrp (1) = ierror
00219          ierrp (2) = (grid_valid_shape(2,1)-grid_valid_shape(1,1)+1) * &
00220                      (grid_valid_shape(2,2)-grid_valid_shape(1,2)+1)
00221 
00222          ierror = PRISM_Error_Alloc
00223          call psmile_error ( ierror, 'vals', ierrp, 2, &
00224                  __FILE__, __LINE__ )
00225          return
00226       endif
00227 !
00228       if (mask_available) then
00229          Allocate (kmask (grid_valid_shape(1,3):grid_valid_shape(2,3)), &
00230                    STAT = ierror)
00231 
00232          if (ierror > 0) then
00233             ierrp (1) = ierror
00234             ierrp (2) = grid_valid_shape(2,3) - grid_valid_shape(1,3) + 1
00235 
00236             ierror = PRISM_Error_Alloc
00237             call psmile_error ( ierror, 'kmask', ierrp, 2, &
00238                  __FILE__, __LINE__ )
00239             return
00240          endif
00241 
00242             do k = grid_valid_shape(1,3), grid_valid_shape(2,3)
00243             kmask (k) = ANY (                                              &
00244                mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1),    &
00245                            grid_valid_shape(1,2):grid_valid_shape(2,2), k) )
00246             end do
00247 
00248 #ifdef PRISM_ASSERTION
00249             if (.not. ANY(kmask)) then
00250                print *, 'grid_valid_shape', grid_valid_shape
00251                call psmile_assert (__FILE__, __LINE__, &
00252                        'at least one entry of kmask should be true')
00253             endif
00254 #endif
00255       endif
00256 !
00257 !===> Brute force
00258 !     TODO: Multigrid search
00259 !
00260       real_huge = huge (distance(1))
00261       k1 = grid_valid_shape (1,3) - 1
00262 !
00263          do n = 1, n_send
00264          if (distance (n) == real_huge) then
00265             dist2 = real_huge
00266          else
00267             dist2 = distance (n) * distance (n)
00268          endif
00269 
00270          fac = real_earth + z_search (n)
00271 !
00272 !===>    ... Compute minimum distance in z-direction
00273 !
00274          if (mask_available) then
00275             kmin = minloc (abs(z_coords (grid_valid_shape(1,3):                  &
00276                                          grid_valid_shape(2,3))) - z_search (n), &
00277                            kmask)
00278          else
00279             kmin = minloc (abs(z_coords (grid_valid_shape(1,3):                  &
00280                                          grid_valid_shape(2,3))) - z_search (n))
00281          endif
00282 !
00283          kmin (1) = kmin (1) + k1
00284 !
00285          distk = abs (z_coords (kmin(1)) - z_search (n))
00286          if (distk >= distance (n)) cycle
00287 !
00288 !===>    ... Compute minimum distance in (lon, lat) direction
00289 ! ???Ist das richtig ???
00290 ! ??? LOOPS auseinandernehmen
00291 
00292 #if 0
00293          do j = grid_valid_shape(1,2), grid_valid_shape(2,2)
00294          vals (:,j) = ( sin_values_lat(j) * sin_search(n, lat)       &
00295                       + cos_values_lat(j) * cos_search(n, lat)       &
00296                         * (cos_values_lon(:) * cos_search(n, lon)    &
00297                           +sin_values_lon(:) * sin_search(n, lon)) )
00298          end do 
00299 #else
00300 !
00301          jend = grid_valid_shape(2,2)
00302 !
00303          vals (:, jend) = (cos_values_lon(:) * cos_search(n, lon)    &
00304                           +sin_values_lon(:) * sin_search(n, lon))
00305 !
00306 !
00307 !===> ...... Note: Loop is split up in 2 part in order
00308 !                  to allow optimization
00309 !
00310 !cdir vector
00311             do j = grid_valid_shape(1,2), jend - 1
00312             vals (:,j) = ( sin_values_lat(j) * sin_search(n, lat)       &
00313                          + cos_values_lat(j) * cos_search(n, lat)       &
00314                            * vals (:,jend))
00315             end do
00316 !
00317 !        j = jend
00318 !
00319          vals (:,j) = ( sin_values_lat(j) * sin_search(n, lat)       &
00320                       + cos_values_lat(j) * cos_search(n, lat)       &
00321                         * vals (:,j))
00322 #endif
00323 !
00324 #ifdef PRISM_ASSERTION
00325          if (minval (vals) < acosm1 - eps .or. &
00326              maxval (vals) > acosp1 + eps) then
00327             print *, 'min vals', minloc (vals), minval (vals)
00328             print *, 'max vals', maxloc (vals), maxval (vals)
00329             call psmile_assert (__FILE__, __LINE__, &
00330                                 'invalid value for acos')
00331          endif
00332 #endif
00333 
00334          vals = min (acosp1, vals)
00335          vals = max (acosm1, vals)
00336 
00337          vals = fac * acos (vals)
00338          vals = vals * vals + distk * distk
00339 !
00340          if (mask_available) then
00341             dist2_mask = dist2
00342             dist_mask  = distance (n)
00343 !
00344 !===> Look for an valid entry
00345 !
00346             if (kmask (kmin(1))) then
00347                ijmin = minloc (vals,                                            &
00348                        mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00349                                    grid_valid_shape(1,2):grid_valid_shape(2,2), &
00350                                    kmin (1)) )
00351                ijmin = ijmin + grid_valid_shape (1, 1:lat) - 1
00352                if (vals (ijmin(1), ijmin(2)) < dist2_mask) then
00353                   dist2_mask = vals (ijmin(1), ijmin(2))
00354                   dist_mask  = sqrt (dist_mask)
00355                endif
00356 !
00357                kex = kmin (1) ! exclude kmin(1); already computed
00358                kmask (kex) = .false.
00359             else
00360                kex = grid_valid_shape(1,3) - 1
00361             endif
00362 !
00363                do k = grid_valid_shape(1,3), grid_valid_shape(2,3)
00364                if (.not. kmask(k)) cycle
00365                if (abs (z_coords (k) - z_search (n)) >= dist_mask) cycle
00366 
00367                   ij = minloc (vals,                                          &
00368                      mask_array (grid_valid_shape(1,1):grid_valid_shape(2,1), &
00369                                  grid_valid_shape(1,2):grid_valid_shape(2,2), &
00370                                  k) )
00371                   ij = ij + grid_valid_shape (1, 1:lat) - 1
00372 
00373                   if (vals (ij(1),ij(2)) < dist_mask) then
00374                      dist2_mask = vals (ijmin(1),ijmin(2))
00375                      dist_mask  = sqrt (dist2_mask)
00376                      ijmin = ij
00377                      kmin (1) = k
00378                   endif
00379                end do
00380 !
00381             if (kex /= grid_valid_shape(1,3)-1) kmask (kex) = .true.
00382 
00383          else
00384 
00385             ijmin = minloc (vals)
00386             ijmin = ijmin + grid_valid_shape (1, 1:lat) - 1
00387 
00388          endif
00389 
00390 !
00391          if (vals(ijmin(1), ijmin(2)) >= dist2) cycle
00392 !
00393 !===>    ... Point with smaller distance found
00394 !
00395          locations (1, n) = ijmin (1)
00396          locations (2, n) = ijmin (2)
00397          locations (3, n) = kmin (1)
00398 !
00399          nfound (n) = 1
00400 !
00401 #ifdef DEBUGX
00402          print *, 'psmile_search_nn_3d_reg_real: n, old, dist2', &
00403                   n, dist2, vals (ijmin(1), ijmin(2))
00404          print *, 'psmile_search_nn_3d_reg_real: n, loc', &
00405                   n, locations (:, n)
00406 #endif
00407 
00408          distance (n) = sqrt (vals (ijmin(1), ijmin(2)))
00409             
00410          enddo ! n
00411 !
00412       Deallocate (vals)
00413 !
00414       if (mask_available) Deallocate (kmask)
00415 !
00416 !===> All done
00417 !
00418 #ifdef VERBOSE
00419       print 9980, trim(ch_id), ierror
00420 
00421       call psmile_flushstd
00422 #endif /* VERBOSE */
00423 !
00424       return
00425 !
00426 !  Formats:
00427 !
00428 #ifdef VERBOSE
00429 9990 format (1x, a, ': psmile_search_nn_3d_reg_real: n_send', i7)
00430 9980 format (1x, a, ': psmile_search_nn_3d_reg_real: eof  ierror =', i4)
00431 #endif /* VERBOSE */
00432 
00433       end subroutine PSMILe_Search_nn_3d_reg_real