psmile_mg_final_gauss2_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, NEC Europe Ltd., London, UK.
00003 ! Copyright 2011, DKRZ, Hamburg, Germany.
00004 ! All rights reserved. Use is subject to OASIS4 license terms.
00005 !-----------------------------------------------------------------------
00006 !BOP
00007 !
00008 ! !ROUTINE: psmile_mg_final_gauss2_dble
00009 !
00010 ! !INTERFACE:
00011 
00012 subroutine psmile_mg_final_gauss2_dble (                                       &
00013                   grid_id, found, locations, fnd_loc_range,                    &
00014                   tgt_coords_x, tgt_coords_y, tgt_coords_shape, search_range,  &
00015                   src_corners_x, src_corners_y, src_corner_shape, nbr_corners, &
00016                   ierror)
00017 
00018 !
00019 ! !USES:
00020 !
00021    use prism_constants
00022    use psmile, dummy_interface => psmile_mg_final_gauss2_dble
00023    use psmile_grid_reduced_gauss
00024 #ifdef DEBUG_TRACE
00025    use psmile_debug_trace
00026 #endif
00027 
00028    implicit none
00029 !
00030 ! !INPUT PARAMETERS:
00031 
00032    integer, intent (in)            :: grid_id
00033 
00034 !     Grid id
00035 
00036    integer, intent (in)            :: fnd_loc_range (2, ndim_3d)
00037 
00038 !     Dimension of "loc" and "found"
00039 
00040    integer, intent (in)            :: tgt_coords_shape (2, ndim_3d)
00041 
00042 !     Dimension of coordinate arrays "coord1", "coord2"
00043 !     which contain the coordinates to be searched.
00044 
00045    double precision, intent (in)   ::                            
00046       tgt_coords_x (tgt_coords_shape(1,1):tgt_coords_shape(2,1), 
00047                     tgt_coords_shape(1,2):tgt_coords_shape(2,2), 
00048                     tgt_coords_shape(1,3):tgt_coords_shape(2,3))
00049 
00050    double precision, intent (in)   ::                            
00051       tgt_coords_y (tgt_coords_shape(1,1):tgt_coords_shape(2,1), 
00052                     tgt_coords_shape(1,2):tgt_coords_shape(2,2), 
00053                     tgt_coords_shape(1,3):tgt_coords_shape(2,3))
00054 
00055 !     Coordinates to be searched
00056 
00057    integer, intent (in)            :: search_range (2, ndim_3d)
00058 
00059 !     Index range in "src_coords" to be searched
00060 
00061    integer, intent (in)            :: src_corner_shape (2, ndim_2d)
00062 
00063 !     Dimension of corner coordinate arrays
00064 
00065    integer, intent (in)            :: nbr_corners
00066 
00067 !     Number of corners
00068 !     Note: This routine assumes nbr_corners == 2.
00069 
00070    double precision,   intent (in) ::                 
00071       src_corners_x ( src_corner_shape(1,1):src_corner_shape(2,1), nbr_corners)
00072 
00073 !     Coordinates of corner points of control volume (Longitude)
00074 
00075    double precision,   intent (in) ::                 
00076       src_corners_y ( src_corner_shape(1,1):src_corner_shape(2,1), nbr_corners)
00077 
00078 !     Coordinates of corner points of control volume (Latitude)
00079 
00080 
00081    integer, intent (inout)         :: found (fnd_loc_range(1,1):fnd_loc_range(2,1), 
00082                                              fnd_loc_range(1,2):fnd_loc_range(2,2), 
00083                                              fnd_loc_range(1,3):fnd_loc_range(2,3))
00084 
00085 !     Finest level number on which a grid cell was found for point i,j,k.
00086 !     Input :
00087 !     Level number < -1: Point was not found and
00088 !                        and last level number was (-found(i,j,k))
00089 !     Level number = -(nlev+1): Never found
00090 !     Output :
00091 !     found(i,j,k) = ??:
00092 !     Otherwise, not contained in the corner grid.
00093 
00094    integer, intent (inout)         :: locations (fnd_loc_range(1,1):fnd_loc_range(2,1), 
00095                                                  fnd_loc_range(1,2):fnd_loc_range(2,2), 
00096                                                  fnd_loc_range(1,3):fnd_loc_range(2,3))
00097 !
00098 !     Indices of the grid cell in which the point was found.
00099 !     The indices are relative to corner shape "src_corner_shape" !
00100 !     I.e. the indices are already transformed into corner/method dimensions.
00101 !     Note: The indices are indices in the Gauss Reduced Grid.
00102 !
00103 ! !OUTPUT PARAMETERS:
00104 !
00105    integer, intent (out)           :: ierror
00106 !
00107 ! !LOCAL VARIABLES
00108 !
00109 !  some loop variables
00110 !
00111    integer                         :: i, j, k, n
00112 !
00113 !  contains the coordinates of a target point
00114 !
00115    type (point_dble)               :: tgt_point
00116 !
00117 !  contains a point and its 8 direct neighbours (1d local indices)
00118 !     7   8   9
00119 !
00120 !     4   5   6
00121 !
00122 !     1   2   3
00123 !
00124    integer                         :: neighbours (9)
00125 !
00126 #ifdef DEBUG_TRACE
00127    logical                         :: debug_trace_on
00128 #endif
00129 !
00130 ! !DEFINED PARAMETERS:
00131 !
00132    integer, parameter              :: not_found = -2
00133 !
00134 ! !DESCRIPTION:
00135 !
00136 ! Subroutine "psmile_mg_final_gauss2_dble" at this point the locations found
00137 ! are just an approximation (due to use of auxiliary grid and remapping to
00138 ! actual reduced gauss grid). This routine does the final step of the search
00139 ! and tries to find the actual cells that match the target coordinates, which
00140 ! were received from the target process.
00141 !
00142 ! !REVISION HISTORY:
00143 !
00144 !   Date      Programmer   Description
00145 ! ----------  ----------   -----------
00146 ! 09.02.11    M. Hanke     created in order to replace old version
00147 !                          of this routine
00148 !
00149 !EOP
00150 !----------------------------------------------------------------------
00151 !
00152 !  $Id: psmile_mg_final_gauss2_dble.F90 2966 2011-02-18 09:47:30Z hanke $
00153 !  $Author: hanke $
00154 !
00155    Character(len=len_cvs_string), save :: mycvs = 
00156        '$Id: psmile_mg_final_gauss2_dble.F90 2966 2011-02-18 09:47:30Z hanke $'
00157 !
00158 !----------------------------------------------------------------------
00159 !
00160 !  Initialization
00161 !
00162    ierror = 0
00163 
00164 
00165 #ifdef VERBOSE
00166    print 9990, trim(ch_id), search_range
00167 
00168    call psmile_flushstd
00169 #endif /* VERBOSE */
00170 
00171 #ifdef PRISM_ASSERTION
00172    if (search_range (1,1) < fnd_loc_range(1,1) .or. &
00173        search_range (2,1) > fnd_loc_range(2,1) .or. &
00174        search_range (1,2) < fnd_loc_range(1,2) .or. &
00175        search_range (2,2) > fnd_loc_range(2,2) .or. &
00176        search_range (1,3) < fnd_loc_range(1,3) .or. &
00177        search_range (2,3) > fnd_loc_range(2,3)) then
00178       call psmile_assert (__FILE__, __LINE__, &
00179                            "search_range exeeds the extent of fnd_loc_range")
00180    endif
00181 #endif
00182 
00183 
00184    ! for all target points that need to be searched for in the local data
00185    do i = search_range (1,1), search_range (2,1)
00186       do j = search_range (1,2), search_range (2,2)
00187          do k = search_range (1,3), search_range (2,3)
00188 
00189 #ifdef DEBUG_TRACE
00190             if (ictl_ind(1) == i .and. ictl_ind(2) == j) then
00191                print 9970, trim(ch_id), "found", found(i,j,k)
00192                print 9970, trim(ch_id), "location",  locations(i,j,k)
00193                debug_trace_on = .true.
00194             else
00195                debug_trace_on = .false.
00196             endif
00197 #endif
00198             ! if the respective point was not found
00199             if (found(i,j,k) /= 1) cycle
00200 
00201             tgt_point%x = tgt_coords_x(i,j,k)
00202             tgt_point%y = tgt_coords_y(i,j,k)
00203 
00204 #ifdef DEBUG_TRACE
00205             if (ictl_ind(1) == i .and. ictl_ind(2) == j) then
00206                print 9971, trim(ch_id), tgt_point
00207             endif
00208 #endif
00209 
00210             ! check the actual cell that is associated with the
00211             ! location found by the multigrid search (5 )
00212             ! if this cell matches tgt point, then found and location are
00213             ! already correct and we do not need to do anything else of
00214             ! this tgt point
00215             if (check_cell (locations(i,j,k), tgt_point)) cycle
00216 
00217             ! if the target point is not in this cell check all
00218             ! 8 neighbours of the cell
00219             ! get the data for the 8 neighbours from the neighbour lookup table
00220             neighbours (1:3) = grids(grid_id)%reduced_gauss_data%local_1d_stencil_lookup(1:3,locations(i,j,k))
00221             neighbours (4:6) = grids(grid_id)%reduced_gauss_data%local_1d_stencil_lookup(5:7,locations(i,j,k))
00222             neighbours (7:9) = grids(grid_id)%reduced_gauss_data%local_1d_stencil_lookup(9:11,locations(i,j,k))
00223             where ( neighbours(:) < grids(grid_id)%grid_shape(1,1) .or. &
00224                     neighbours(:) > grids(grid_id)%grid_shape(2,1))
00225 
00226                neighbours(:) = psmile_undef
00227             endwhere
00228 
00229             do n = 1, 9
00230 
00231                ! we already checked n == 5 (== locations(i,j,k))
00232                if (n == 5) cycle
00233 
00234                ! neighbours(n) == psmile_undef means that the cell is not in our local partition
00235                if (neighbours(n) == psmile_undef) cycle
00236 
00237                ! check whether the current neighour cell contains the target point
00238                if (check_cell (neighbours(n), tgt_point)) then
00239 
00240                   ! we need to adjust the location array because the point was
00241                   ! found in a neighbour of the cell the current location is
00242                   ! pointing to
00243                   locations (i,j,k) = neighbours(n)
00244 
00245                   exit
00246                endif
00247 
00248             enddo ! n
00249 
00250             ! if no matching source cell was found, we set the target point to not_found
00251             if (n > 9) found (i,j,k) = not_found
00252 
00253 #ifdef DEBUG_TRACE
00254             if (ictl_ind(1) == i .and. ictl_ind(2) == j) then
00255                print 9970, trim(ch_id), "new found", found (i,j,k)
00256                print 9970, trim(ch_id), "new location", locations (i,j,k)
00257             endif
00258 #endif
00259          enddo ! k
00260       enddo ! j
00261    enddo ! i
00262 !
00263 !===> All done
00264 !
00265 !
00266 #ifdef VERBOSE
00267    print 9980, trim(ch_id)
00268 
00269    call psmile_flushstd
00270 #endif /* VERBOSE */
00271 !
00272 !  Formats:
00273 !
00274 9990  format (1x, a, ': psmile_mg_final_gauss2_dble:  search_range ', 6i10)
00275 9980  format (1x, a, ': psmile_mg_final_gauss2_dble: eof')
00276 
00277 #ifdef DEBUG_TRACE
00278 9970  format (1x, a, ': psmile_mg_final_gauss2_dble: ictl ', a, ':', 1i10)
00279 9971  format (1x, a, ': psmile_mg_final_gauss2_dble: ictl tgt_point:', 2e13.6)
00280 #endif
00281 
00282    contains
00283 
00284       logical function check_cell (src_cell_loc, tgt_point)
00285 
00286          use psmile_grid, only : common_grid_range
00287 
00288          integer, intent (in)           :: src_cell_loc
00289          type (point_dble), intent (in) :: tgt_point
00290 
00291          type (point_dble)              :: tgt_point_
00292 
00293          tgt_point_ = tgt_point
00294 
00295          do while (tgt_point_%x < minval(src_corners_x (src_cell_loc, 1:2)))
00296             tgt_point_%x = tgt_point_%x + common_grid_range(2,1) - common_grid_range(1,1)
00297          enddo
00298 
00299          do while (tgt_point_%x > maxval(src_corners_x (src_cell_loc, 1:2)))
00300             tgt_point_%x = tgt_point_%x - (common_grid_range(2,1) - common_grid_range(1,1))
00301          enddo
00302 
00303          check_cell = tgt_point_%x <= maxval(src_corners_x (src_cell_loc, 1:2)) .and. &
00304                       tgt_point_%x >= minval(src_corners_x (src_cell_loc, 1:2)) .and. &
00305                       tgt_point_%y <= maxval(src_corners_y (src_cell_loc, 1:2)) .and. &
00306                       tgt_point_%y >= minval(src_corners_y (src_cell_loc, 1:2))
00307 
00308 #ifdef DEBUG_TRACE
00309          if (debug_trace_on) then
00310             print 9972, trim(ch_id), "testing source cell:"
00311             print "(2e13.6)", src_corners_x (src_cell_loc, 1:2)
00312             print "(2e13.6)", src_corners_y (src_cell_loc, 1:2)
00313             if (check_cell) then
00314                print 9972, trim(ch_id), "point is in this cell"
00315             else
00316                print 9972, trim(ch_id), "point is not in this cell"
00317             endif
00318          endif
00319 #endif
00320 9972  format (1x, a, ': psmile_mg_final_gauss2_dble: ictl ', a)
00321       end function check_cell
00322 
00323 end subroutine psmile_mg_final_gauss2_dble