geogrid/src/proc_point_module.F

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! Module: proc_point_module
   3 !
   4 ! Purpose: This module provides routines that produce a value for a model grid
   5 !    point in two ways. If the field for which a value is being calculated is
   6 !    a continuous field, this module provided functionality to interpolate
   7 !    from the source array to the specified point. If the field is a categorical
   8 !    field, this module provided functionality to accumulate the values of all
   9 !    source points whose nearest model gridpoint is the specified point.
  10 !    Routines are also provided that help the caller determine an optimized
  11 !    order in which to process the model grid points.
  12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  13 module proc_point_module
  14
  15    ! Modules
  16    use bitarray_module
  17    use hash_module
  18    use misc_definitions_module
  19    use module_debug
  20    use source_data_module
  21
  22    ! Information about which tile is in memory
  23    integer :: src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, src_npts_bdr
  24    integer :: src_level
  25    character (len=128) :: src_fieldname
  26    character (len=256) :: src_fname
  27
  28    ! Source tiles
  29    real, pointer, dimension(:,:,:) :: src_array
  30
  31    ! Hash to track which tiles we have already processed
  32    type (hashtable) :: h_table
  33
  34    contains
  35
  36    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  37    ! Name: proc_point_init
  38    !
  39    ! Purpose: Initialize the module.
  40    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  41    subroutine proc_point_init()
  42
  43       implicit none
  44
  45       ! Initialize module variables
  46       src_min_x = INVALID
  47       src_min_y = INVALID
  48       src_min_z = INVALID
  49       src_max_x = INVALID
  50       src_max_y = INVALID
  51       src_max_z = INVALID
  52       src_fieldname = ' '
  53       src_fname = ' '
  54       nullify(src_array)
  55
  56       call hash_init(h_table)
  57
  58    end subroutine proc_point_init
  59
  60
  61    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  62    ! Name: proc_point_shutdown
  63    !
  64    ! Purpose: Do any cleanup work.
  65    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  66    subroutine proc_point_shutdown()
  67
  68       implicit none
  69
  70       ! Effectively reset the hash table that tracks which tiles have been processed
  71       !   by removing all entries
  72       call hash_destroy(h_table)
  73
  74       if (associated(src_array)) deallocate(src_array)
  75
  76       src_min_x = INVALID
  77       src_min_y = INVALID
  78       src_min_z = INVALID
  79       src_max_x = INVALID
  80       src_max_y = INVALID
  81       src_max_z = INVALID
  82       src_fieldname = ' '
  83       src_fname = ' '
  84
  85    end subroutine proc_point_shutdown
  86
  87
  88    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  89    ! Name: accum_categorical
  90    !
  91    ! Purpose: Count the number of source points in each category whose nearest
  92    !   neighbor is the specified model grid point.
  93    !
  94    ! NOTE: When processing the source tile, those source points that are
  95    !   closest to a different model grid point will be added to the totals for
  96    !   such grid points; thus, an entire source tile will be processed at a time.
  97    !   This routine really processes for all model grid points that are
  98    !   within a source tile, and not just for a single grid point.
  99    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 100    subroutine accum_categorical(xlat, xlon, istagger, array, &
 101                                 start_i, end_i, start_j, end_j, &
 102                                 start_k, end_k, fieldname, processed_pts, &
 103                                 new_pts, ilevel, msgval, maskval, sr_x, sr_y)
 104
 105       use llxy_module
 106       use bitarray_module
 107
 108       implicit none
 109
 110       ! Arguments
 111       integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, &
 112                              istagger, ilevel
 113       real, intent(in) :: xlat, xlon, msgval, maskval
 114       real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: array
 115       character (len=128), intent(in) :: fieldname
 116       type (bitarray), intent(inout) :: processed_pts, new_pts
 117       integer, intent(in), optional :: sr_x, sr_y
 118
 119       ! Local variables
 120       integer :: istatus, i, j
 121       integer :: current_domain, k
 122       integer, pointer, dimension(:,:,:) :: where_maps_to
 123       real :: rlat, rlon
 124       real :: rarea
 125       real :: rsr_x, rsr_y
 126
 127       rlat = xlat
 128       if (xlon >= 180.) then
 129          rlon = xlon - 360.
 130       else
 131          rlon = xlon
 132       end if
 133
 134       rsr_x = 1.0
 135       rsr_y = 1.0
 136       if (present(sr_x)) rsr_x = real(sr_x)
 137       if (present(sr_y)) rsr_y = real(sr_y)
 138
 139       ! Assume source data is on unstaggered grid; specify M for istagger argument
 140       call get_data_tile(rlat, rlon, ilevel, fieldname, &
 141                          src_fname, src_array, src_min_x, src_max_x, src_min_y, &
 142                          src_max_y, src_min_z, src_max_z, src_npts_bdr, &
 143                          istatus)
 144
 145       src_fieldname = fieldname
 146       src_level = ilevel
 147
 148       call hash_insert(h_table, src_fname)
 149
 150       if (istatus /= 0) return
 151
 152       allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
 153       do i=src_min_x,src_max_x
 154          do j=src_min_y,src_max_y
 155             where_maps_to(i,j,1) = NOT_PROCESSED
 156          end do
 157       end do
 158
 159       call process_categorical_block(src_array, istagger, where_maps_to, &
 160                               src_min_x+src_npts_bdr, src_min_y+src_npts_bdr, src_min_z, &
 161                               src_max_x-src_npts_bdr, src_max_y-src_npts_bdr, src_max_z, &
 162                               array, start_i, end_i, start_j, end_j, start_k, end_k, &
 163                               processed_pts, new_pts, ilevel, rsr_x, rsr_y, msgval, maskval)
 164
 165       ! If a grid cell has less than half of its area covered by data from this source,
 166       !   then clear the cell and let another source fill in the cell
 167       if (ilevel > 1) then
 168          do i=start_i,end_i
 169             do j=start_j,end_j
 170                if (bitarray_test(new_pts, i-start_i+1, j-start_j+1) .and. &
 171                    .not. bitarray_test(processed_pts, i-start_i+1, j-start_j+1)) then
 172                   rarea = 0.
 173                   do k=start_k,end_k
 174                      rarea = rarea + array(i,j,k)
 175                   end do
 176                   current_domain = iget_selected_domain()
 177                   call select_domain(SOURCE_PROJ)
 178                   if (proj_stack(current_nest_number)%dx < 0.) then
 179                      rarea = rarea * (proj_stack(current_nest_number)%latinc*111000.)**2.0
 180                   else
 181                      rarea = rarea * proj_stack(current_nest_number)%dx**2.0
 182                   end if
 183                   call select_domain(current_domain)
 184                   if (proj_stack(current_nest_number)%dx < 0.) then
 185                      if ((proj_stack(current_nest_number)%latinc*111000.)**2.0 > 2.0*rarea) then
 186                         do k=start_k,end_k
 187                            array(i,j,k) = 0.
 188                         end do
 189                         call bitarray_clear(new_pts, i-start_i+1, j-start_j+1)
 190                      end if
 191                   else
 192                      if (proj_stack(current_nest_number)%dx**2.0 > 2.0*rarea) then
 193                         do k=start_k,end_k
 194                            array(i,j,k) = 0.
 195                         end do
 196                         call bitarray_clear(new_pts, i-start_i+1, j-start_j+1)
 197                      end if
 198                   end if
 199                end if
 200             end do
 201          end do
 202       end if
 203
 204       deallocate(where_maps_to)
 205
 206    end subroutine accum_categorical
 207
 208
 209    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 210    ! Name: process_categorical_block
 211    !
 212    ! Purpose: To recursively process a subarray of categorical data, assigning
 213    !   the points in a block to their nearest grid point. The nearest neighbor
 214    !   may be estimated in some cases; for example, if the four corners of a
 215    !   subarray all have the same nearest grid point, all elements in the
 216    !   subarray are assigned to that grid point.
 217    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 218    recursive subroutine process_categorical_block(tile_array, istagger, where_maps_to, &
 219                                    min_i, min_j, min_k, max_i, max_j, max_k, dst_array, &
 220                                    start_x, end_x, start_y, end_y, start_z, end_z, &
 221                                    processed_pts, new_pts, ilevel, sr_x, sr_y, &
 222                                    msgval, maskval, mask_array)
 223
 224       use llxy_module
 225
 226       implicit none
 227
 228       ! Arguments
 229       integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, istagger, &
 230                              start_x, end_x, start_y, end_y, start_z, end_z, ilevel
 231       integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
 232       real, intent(in) :: sr_x, sr_y, msgval, maskval
 233       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
 234       real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array
 235       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
 236       type (bitarray), intent(inout) :: processed_pts, new_pts
 237
 238       ! Local variables
 239       integer :: x_dest, y_dest, i, j, k, center_i, center_j, current_domain
 240       real :: lat_corner, lon_corner, rx, ry
 241
 242       ! Compute the model grid point that the corners of the rectangle to be
 243       !   processed map to
 244       ! Lower-left corner
 245       if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
 246          current_domain = iget_selected_domain()
 247          call select_domain(SOURCE_PROJ)
 248          call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, M)
 249          call select_domain(current_domain)
 250          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 251          rx = (rx-1.0) * sr_x + 1.0
 252          ry = (ry-1.0) * sr_y + 1.0
 253          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 254             where_maps_to(min_i,min_j,1) = nint(rx)
 255             where_maps_to(min_i,min_j,2) = nint(ry)
 256          else
 257             where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
 258          end if
 259       end if
 260
 261       ! Upper-left corner
 262       if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
 263          current_domain = iget_selected_domain()
 264          call select_domain(SOURCE_PROJ)
 265          call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, M)
 266          call select_domain(current_domain)
 267          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 268          rx = (rx-1.0) * sr_x + 1.0
 269          ry = (ry-1.0) * sr_y + 1.0
 270          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 271             where_maps_to(min_i,max_j,1) = nint(rx)
 272             where_maps_to(min_i,max_j,2) = nint(ry)
 273          else
 274             where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
 275          end if
 276       end if
 277
 278       ! Upper-right corner
 279       if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
 280          current_domain = iget_selected_domain()
 281          call select_domain(SOURCE_PROJ)
 282          call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, M)
 283          call select_domain(current_domain)
 284          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 285          rx = (rx-1.0) * sr_x + 1.0
 286          ry = (ry-1.0) * sr_y + 1.0
 287          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 288             where_maps_to(max_i,max_j,1) = nint(rx)
 289             where_maps_to(max_i,max_j,2) = nint(ry)
 290          else
 291             where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
 292          end if
 293       end if
 294
 295       ! Lower-right corner
 296       if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
 297          current_domain = iget_selected_domain()
 298          call select_domain(SOURCE_PROJ)
 299          call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, M)
 300          call select_domain(current_domain)
 301          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 302          rx = (rx-1.0) * sr_x + 1.0
 303          ry = (ry-1.0) * sr_y + 1.0
 304          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 305             where_maps_to(max_i,min_j,1) = nint(rx)
 306             where_maps_to(max_i,min_j,2) = nint(ry)
 307          else
 308             where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
 309          end if
 310       end if
 311
 312       ! If all four corners map to same model grid point, accumulate the
 313       !   entire rectangle
 314       if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
 315           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
 316           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
 317           where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
 318           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
 319           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
 320           where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
 321          x_dest = where_maps_to(min_i,min_j,1)
 322          y_dest = where_maps_to(min_i,min_j,2)
 323
 324          ! If this grid point was already given a value from higher-priority source data,
 325          !   there is nothing to do.
 326          if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 327
 328             ! If this grid point has never been given a value by this level of source data,
 329             !   initialize the point
 330             if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 331                do k=start_z,end_z
 332                   dst_array(x_dest,y_dest,k) = 0.
 333                end do
 334             end if
 335
 336             ! Count all the points whose nearest neighbor is this grid point
 337             if (present(mask_array)) then
 338                do i=min_i,max_i
 339                   do j=min_j,max_j
 340                      ! Ignore masked/missing values in the source data
 341                      if ((tile_array(i,j,min_k) /= msgval) .and. &
 342                          (mask_array(i,j) /= maskval)) then
 343                         if (int(tile_array(i,j,min_k)) >= start_z .and. int(tile_array(i,j,min_k)) <= end_z) then
 344                            dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) = &
 345                                dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) + 1.0
 346                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 347                         else
 348                            call mprintf(.true., WARN, 'In source tile %s, point (%i, %i) has '// &
 349                                         'an invalid category of %i', &
 350                                         s1=trim(src_fname), i1=i, i2=j, i3=int(tile_array(i,j,min_k)))
 351                         end if
 352                      end if
 353                   end do
 354                end do
 355             else
 356                do i=min_i,max_i
 357                   do j=min_j,max_j
 358                      ! Ignore masked/missing values in the source data
 359                      if (tile_array(i,j,min_k) /= msgval) then
 360                         if (int(tile_array(i,j,min_k)) >= start_z .and. int(tile_array(i,j,min_k)) <= end_z) then
 361                            dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) = &
 362                                dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) + 1.0
 363                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 364                         else
 365                            call mprintf(.true., WARN, 'In source tile %s, point (%i, %i) '// &
 366                                         'has an invalid category of %i', &
 367                                         s1=trim(src_fname), i1=i, i2=j, i3=int(tile_array(i,j,min_k)))
 368                         end if
 369                      end if
 370                   end do
 371                end do
 372             end if
 373
 374          end if
 375
 376       ! Rectangle is a square of four points, and we can simply deal with each of the points
 377       else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
 378          do i=min_i,max_i
 379             do j=min_j,max_j
 380                x_dest = where_maps_to(i,j,1)
 381                y_dest = where_maps_to(i,j,2)
 382
 383                if (x_dest /= OUTSIDE_DOMAIN) then
 384
 385                   if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 386                      if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 387                         do k=start_z,end_z
 388                            dst_array(x_dest,y_dest,k) = 0.
 389                         end do
 390                      end if
 391
 392                      ! Ignore masked/missing values
 393                      if (present(mask_array)) then
 394                         if ((tile_array(i,j,min_k) /= msgval) .and. &
 395                              (mask_array(i,j) /= maskval)) then
 396                            if (int(tile_array(i,j,min_k)) >= start_z .and. int(tile_array(i,j,min_k)) <= end_z) then
 397                               dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) = &
 398                                   dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) + 1.0
 399                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 400                            else
 401                               call mprintf(.true., WARN, 'In source tile %s, point (%i, %i) has '// &
 402                                            'an invalid category of %i', &
 403                                            s1=trim(src_fname), i1=i, i2=j, i3=int(tile_array(i,j,min_k)))
 404                            end if
 405                         end if
 406                      else
 407                         if (tile_array(i,j,min_k) /= msgval) then
 408                            if (int(tile_array(i,j,min_k)) >= start_z .and. int(tile_array(i,j,min_k)) <= end_z) then
 409                               dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) = &
 410                                   dst_array(x_dest,y_dest,int(tile_array(i,j,min_k))) + 1.0
 411                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 412                            else
 413                               call mprintf(.true., WARN, 'In source tile %s, point (%i, %i) has '// &
 414                                            'an invalid category of %i', &
 415                                            s1=trim(src_fname), i1=i, i2=j, i3=int(tile_array(i,j,min_k)))
 416                            end if
 417                         end if
 418                      end if
 419                   end if
 420
 421                end if
 422             end do
 423          end do
 424
 425       ! Not all corners map to the same grid point, and the rectangle contains more than
 426       !   four points
 427       else
 428          center_i = (max_i + min_i)/2
 429          center_j = (max_j + min_j)/2
 430
 431          ! Recursively process lower-left rectangle
 432          call process_categorical_block(tile_array, istagger, where_maps_to, min_i, min_j, min_k, center_i, &
 433                     center_j, max_k, dst_array, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 434                     new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 435
 436          ! Recursively process lower-right rectangle
 437          if (center_i < max_i) then
 438             call process_categorical_block(tile_array, istagger, where_maps_to, center_i+1, min_j, min_k, max_i, &
 439                        center_j, max_k, dst_array, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 440                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 441          end if
 442
 443          ! Recursively process upper-left rectangle
 444          if (center_j < max_j) then
 445             call process_categorical_block(tile_array, istagger, where_maps_to, min_i, center_j+1, min_k, center_i, &
 446                        max_j, max_k, dst_array, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 447                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 448          end if
 449
 450          ! Recursively process upper-right rectangle
 451          if (center_i < max_i .and. center_j < max_j) then
 452             call process_categorical_block(tile_array, istagger, where_maps_to, center_i+1, center_j+1, min_k, max_i, &
 453                        max_j, max_k, dst_array, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 454                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 455          end if
 456       end if
 457
 458    end subroutine process_categorical_block
 459
 460
 461    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 462    ! Name: accum_continuous
 463    !
 464    ! Purpose: Sum up all of the source data points whose nearest neighbor in the
 465    !   model grid is the specified model grid point.
 466    !
 467    ! NOTE: When processing the source tile, those source points that are
 468    !   closest to a different model grid point will be added to the totals for
 469    !   such grid points; thus, an entire source tile will be processed at a time.
 470    !   This routine really processes for all model grid points that are
 471    !   within a source tile, and not just for a single grid point.
 472    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 473    subroutine accum_continuous(xlat, xlon, istagger, array, n, &
 474                                 start_i, end_i, start_j, end_j, &
 475                                 start_k, end_k, fieldname, processed_pts, &
 476                                 new_pts, ilevel, msgval, maskval, sr_x, sr_y)
 477
 478       implicit none
 479
 480       ! Arguments
 481       integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, &
 482                              istagger, ilevel
 483       real, intent(in) :: xlat, xlon, msgval, maskval
 484       real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: array, n
 485       character (len=128), intent(in) :: fieldname
 486       type (bitarray), intent(inout) :: processed_pts, new_pts
 487       integer, intent(in), optional :: sr_x, sr_y
 488
 489       ! Local variables
 490       integer :: istatus, i, j
 491       integer, pointer, dimension(:,:,:) :: where_maps_to
 492       real :: rlat, rlon
 493       real :: rsr_x, rsr_y
 494
 495       rlat = xlat
 496       if (xlon >= 180.) then
 497          rlon = xlon - 360.
 498       else
 499          rlon = xlon
 500       end if
 501
 502       rsr_x = 1.0
 503       rsr_y = 1.0
 504       if (present(sr_x)) rsr_x = real(sr_x)
 505       if (present(sr_y)) rsr_y = real(sr_y)
 506
 507       ! Assume source data is on unstaggered grid; specify M for istagger argument
 508       call get_data_tile(rlat, rlon, ilevel, fieldname, &
 509                          src_fname, src_array, src_min_x, src_max_x, src_min_y, &
 510                          src_max_y, src_min_z, src_max_z, src_npts_bdr, &
 511                          istatus)
 512
 513       src_fieldname = fieldname
 514       src_level = ilevel
 515
 516       call hash_insert(h_table, src_fname)
 517
 518       if (istatus /= 0) then
 519          src_min_x = INVALID
 520          src_min_y = INVALID
 521          src_max_x = INVALID
 522          src_max_y = INVALID
 523          return
 524       end if
 525
 526       allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
 527       do i=src_min_x,src_max_x
 528          do j=src_min_y,src_max_y
 529             where_maps_to(i,j,1) = NOT_PROCESSED
 530          end do
 531       end do
 532
 533       call process_continuous_block(src_array, istagger, where_maps_to, &
 534                                src_min_x+src_npts_bdr, src_min_y+src_npts_bdr, src_min_z, &
 535                                src_max_x-src_npts_bdr, src_max_y-src_npts_bdr, src_max_z, &
 536                                array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
 537                                processed_pts, new_pts, ilevel, rsr_x, rsr_y, msgval, maskval)
 538
 539       deallocate(where_maps_to)
 540
 541    end subroutine accum_continuous
 542
 543
 544    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 545    ! Name: process_continuous_block
 546    !
 547    ! Purpose: To recursively process a subarray of continuous data, adding the
 548    !   points in a block to the sum for their nearest grid point. The nearest
 549    !   neighbor may be estimated in some cases; for example, if the four corners
 550    !   of a subarray all have the same nearest grid point, all elements in the
 551    !   subarray are added to that grid point.
 552    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 553    recursive subroutine process_continuous_block(tile_array, istagger, where_maps_to, &
 554                                    min_i, min_j, min_k, max_i, max_j, max_k, dst_array, n, &
 555                                    start_x, end_x, start_y, end_y, start_z, end_z, &
 556                                    processed_pts, new_pts, ilevel, sr_x, sr_y, &
 557                                    msgval, maskval, mask_array)
 558
 559       use llxy_module
 560
 561       implicit none
 562
 563       ! Arguments
 564       integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, istagger, &
 565                              start_x, end_x, start_y, end_y, start_z, end_z, ilevel
 566       integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
 567       real, intent(in) :: sr_x, sr_y, msgval, maskval
 568       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
 569       real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
 570       real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
 571       type (bitarray), intent(inout) :: processed_pts, new_pts
 572
 573       ! Local variables
 574       integer :: x_dest, y_dest, i, j, k, center_i, center_j, current_domain
 575       real :: lat_corner, lon_corner, rx, ry
 576
 577       ! Compute the model grid point that the corners of the rectangle to be
 578       !   processed map to
 579       ! Lower-left corner
 580       if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
 581          current_domain = iget_selected_domain()
 582          call select_domain(SOURCE_PROJ)
 583          call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, M)
 584          call select_domain(current_domain)
 585          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 586          rx = (rx-1.0) * sr_x + 1.0
 587          ry = (ry-1.0) * sr_y + 1.0
 588          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 589             where_maps_to(min_i,min_j,1) = nint(rx)
 590             where_maps_to(min_i,min_j,2) = nint(ry)
 591          else
 592             where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
 593          end if
 594       end if
 595
 596       ! Upper-left corner
 597       if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
 598          current_domain = iget_selected_domain()
 599          call select_domain(SOURCE_PROJ)
 600          call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, M)
 601          call select_domain(current_domain)
 602          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 603          rx = (rx-1.0) * sr_x + 1.0
 604          ry = (ry-1.0) * sr_y + 1.0
 605          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 606             where_maps_to(min_i,max_j,1) = nint(rx)
 607             where_maps_to(min_i,max_j,2) = nint(ry)
 608          else
 609             where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
 610          end if
 611       end if
 612
 613       ! Upper-right corner
 614       if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
 615          current_domain = iget_selected_domain()
 616          call select_domain(SOURCE_PROJ)
 617          call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, M)
 618          call select_domain(current_domain)
 619          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 620          rx = (rx-1.0) * sr_x + 1.0
 621          ry = (ry-1.0) * sr_y + 1.0
 622          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 623             where_maps_to(max_i,max_j,1) = nint(rx)
 624             where_maps_to(max_i,max_j,2) = nint(ry)
 625          else
 626             where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
 627          end if
 628       end if
 629
 630       ! Lower-right corner
 631       if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
 632          current_domain = iget_selected_domain()
 633          call select_domain(SOURCE_PROJ)
 634          call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, M)
 635          call select_domain(current_domain)
 636          call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
 637          rx = (rx-1.0) * sr_x + 1.0
 638          ry = (ry-1.0) * sr_y + 1.0
 639          if (real(start_x) <= rx .and. rx <= real(end_x) .and. real(start_y) <= ry .and. ry <= real(end_y)) then
 640             where_maps_to(max_i,min_j,1) = nint(rx)
 641             where_maps_to(max_i,min_j,2) = nint(ry)
 642          else
 643             where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
 644          end if
 645       end if
 646
 647       ! If all four corners map to same model grid point, accumulate the
 648       !   entire rectangle
 649       if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
 650           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
 651           where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
 652           where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
 653           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
 654           where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
 655           where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
 656          x_dest = where_maps_to(min_i,min_j,1)
 657          y_dest = where_maps_to(min_i,min_j,2)
 658
 659          ! If this grid point was already given a value from higher-priority source data,
 660          !   there is nothing to do.
 661          if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 662
 663             ! If this grid point has never been given a value by this level of source data,
 664             !   initialize the point
 665             if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 666                do k=min_k,max_k
 667                   dst_array(x_dest,y_dest,k) = 0.
 668                end do
 669             end if
 670
 671             ! Sum all the points whose nearest neighbor is this grid point
 672             if (present(mask_array)) then
 673                do i=min_i,max_i
 674                   do j=min_j,max_j
 675                      do k=min_k,max_k
 676                         ! Ignore masked/missing values in the source data
 677                         if ((tile_array(i,j,k) /= msgval) .and. &
 678                             (mask_array(i,j) /= maskval)) then
 679                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
 680                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
 681                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 682                         end if
 683                      end do
 684                   end do
 685                end do
 686             else
 687                do i=min_i,max_i
 688                   do j=min_j,max_j
 689                      do k=min_k,max_k
 690                         ! Ignore masked/missing values in the source data
 691                         if (tile_array(i,j,k) /= msgval) then
 692                            dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
 693                            n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
 694                            call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 695                         end if
 696                      end do
 697                   end do
 698                end do
 699             end if
 700
 701          end if
 702
 703       ! Rectangle is a square of four points, and we can simply deal with each of the points
 704       else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
 705          do i=min_i,max_i
 706             do j=min_j,max_j
 707                x_dest = where_maps_to(i,j,1)
 708                y_dest = where_maps_to(i,j,2)
 709
 710                if (x_dest /= OUTSIDE_DOMAIN) then
 711
 712                   if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 713                      if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
 714                         do k=min_k,max_k
 715                            dst_array(x_dest,y_dest,k) = 0.
 716                         end do
 717                      end if
 718
 719                      if (present(mask_array)) then
 720                         do k=min_k,max_k
 721                            ! Ignore masked/missing values
 722                            if ((tile_array(i,j,k) /= msgval) .and. &
 723                                 (mask_array(i,j) /= maskval)) then
 724                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
 725                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
 726                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 727                            end if
 728                         end do
 729                      else
 730                         do k=min_k,max_k
 731                            ! Ignore masked/missing values
 732                            if (tile_array(i,j,k) /= msgval) then
 733                               dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
 734                               n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
 735                               call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
 736                            end if
 737                         end do
 738                      end if
 739                   end if
 740
 741                end if
 742             end do
 743          end do
 744
 745       ! Not all corners map to the same grid point, and the rectangle contains more than
 746       !   four points
 747       else
 748          center_i = (max_i + min_i)/2
 749          center_j = (max_j + min_j)/2
 750
 751          ! Recursively process lower-left rectangle
 752          call process_continuous_block(tile_array, istagger, where_maps_to, min_i, min_j, min_k, center_i, &
 753                     center_j, max_k, dst_array, n, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 754                     new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 755
 756          ! Recursively process lower-right rectangle
 757          if (center_i < max_i) then
 758             call process_continuous_block(tile_array, istagger, where_maps_to, center_i+1, min_j, min_k, max_i, &
 759                        center_j, max_k, dst_array, n, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 760                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 761          end if
 762
 763          ! Recursively process upper-left rectangle
 764          if (center_j < max_j) then
 765             call process_continuous_block(tile_array, istagger, where_maps_to, min_i, center_j+1, min_k, center_i, &
 766                        max_j, max_k, dst_array, n, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 767                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 768          end if
 769
 770          ! Recursively process upper-right rectangle
 771          if (center_i < max_i .and. center_j < max_j) then
 772             call process_continuous_block(tile_array, istagger, where_maps_to, center_i+1, center_j+1, min_k, max_i, &
 773                        max_j, max_k, dst_array, n, start_x, end_x, start_y, end_y, start_z, end_z, processed_pts, &
 774                        new_pts, ilevel, sr_x, sr_y, msgval, maskval, mask_array)
 775          end if
 776
 777       end if
 778
 779    end subroutine process_continuous_block
 780
 781
 782    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 783    ! Name: get_point
 784    !
 785    ! Purpose: For a specified lat/lon and level, return the value of the field
 786    !   interpolated to or nearest the lat/lon.
 787    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 788    function get_point(xlat, xlon, lvl, fieldname, &
 789                       ilevel, interp_type, msgval)
 790
 791       ! Modules
 792       use interp_module
 793       use llxy_module
 794
 795       implicit none
 796
 797       ! Arguments
 798       integer, intent(in) :: lvl, ilevel
 799       real, intent(in) :: xlat, xlon, msgval
 800       character (len=128), intent(in) :: fieldname
 801       integer, dimension(:), intent(in) :: interp_type
 802
 803       ! Return value
 804       real :: get_point
 805
 806       ! Local variables
 807       integer :: istatus, current_domain
 808       real :: rlat, rlon, rx, ry
 809
 810       rlat = xlat
 811       if (xlon >= 180.) then
 812          rlon = xlon - 360.
 813       else
 814          rlon = xlon
 815       end if
 816
 817       ! If tile is in memory, interpolate
 818       if (ilevel == src_level .and. is_point_in_tile(rlat, rlon, ilevel) .and. fieldname == src_fieldname) then
 819
 820          current_domain = iget_selected_domain()
 821          call select_domain(SOURCE_PROJ)
 822          call lltoxy(rlat, rlon, rx, ry, M)  ! Assume source data on unstaggered grid
 823          call select_domain(current_domain)
 824
 825          get_point = interp_sequence(rx, ry, lvl, src_array, src_min_x, src_max_x, src_min_y, &
 826                                      src_max_y, src_min_z, src_max_z, msgval, interp_type, 1)
 827
 828       else
 829
 830          call get_data_tile(rlat, rlon, ilevel, fieldname, &
 831                             src_fname, src_array, src_min_x, src_max_x, src_min_y, &
 832                             src_max_y, src_min_z, src_max_z, src_npts_bdr, &
 833                             istatus)
 834
 835          src_fieldname = fieldname
 836          src_level = ilevel
 837
 838          if (istatus /= 0) then
 839             get_point = msgval
 840             return
 841          end if
 842
 843          current_domain = iget_selected_domain()
 844          call select_domain(SOURCE_PROJ)
 845          call lltoxy(rlat, rlon, rx, ry, M)  ! Assume source data on unstaggered grid
 846          call select_domain(current_domain)
 847
 848          get_point = interp_sequence(rx, ry, lvl, src_array, src_min_x, src_max_x, src_min_y, &
 849                                      src_max_y, src_min_z, src_max_z, msgval, interp_type, 1)
 850       end if
 851
 852       return
 853
 854    end function get_point
 855
 856
 857    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 858    ! Name: have_processed_tile
 859    !
 860    ! Purpose: This funtion returns .true. if the tile of data for
 861    !   the specified field has already been processed, and .false. otherwise.
 862    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 863    function have_processed_tile(xlat, xlon, fieldname, ilevel)
 864
 865       implicit none
 866
 867       ! Arguments
 868       integer, intent(in) :: ilevel
 869       real, intent(in) :: xlat, xlon
 870       character (len=128), intent(in) :: fieldname
 871
 872       ! Return value
 873       logical :: have_processed_tile
 874
 875       ! Local variables
 876       integer :: istatus
 877       character (len=256) :: test_fname
 878
 879       call get_tile_fname(test_fname, xlat, xlon, ilevel, fieldname, istatus)
 880       have_processed_tile = hash_search(h_table, test_fname)
 881
 882       return
 883
 884    end function have_processed_tile
 885
 886
 887    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 888    ! Name: is_point_in_tile
 889    !
 890    ! Purpose: Returns whether the specified lat/lon could be processed
 891    !   without incurring a file access.
 892    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 893    function is_point_in_tile(xlat, xlon, ilevel)
 894
 895       use llxy_module
 896
 897       implicit none
 898
 899       ! Arguments
 900       integer, intent(in) :: ilevel
 901       real, intent(in) :: xlat, xlon
 902
 903       ! Return value
 904       logical :: is_point_in_tile
 905
 906       ! Local variables
 907       integer :: current_domain
 908       real :: rlat, rlon, rx, ry
 909
 910       rlat = xlat
 911       if (xlon >= 180.) then
 912          rlon = xlon - 360.
 913       else
 914          rlon = xlon
 915       end if
 916
 917       current_domain = iget_selected_domain()
 918       call select_domain(SOURCE_PROJ)
 919       call lltoxy(rlat, rlon, rx, ry, M)
 920       call select_domain(current_domain)
 921
 922   !    if (real(src_min_x+src_npts_bdr) <= rx .and. rx <= real(src_max_x-src_npts_bdr) .and. &
 923   !        real(src_min_y+src_npts_bdr) <= ry .and. ry <= real(src_max_y-src_npts_bdr)) then
 924 ! BUG 2006-06-01
 925 !      if (src_min_x+src_npts_bdr <= ceiling(rx) .and. floor(rx) <= src_max_x-src_npts_bdr .and. &
 926 !          src_min_y+src_npts_bdr <= ceiling(ry) .and. floor(ry) <= src_max_y-src_npts_bdr) then
 927       if (src_min_x+src_npts_bdr <= floor(rx+0.5) .and. ceiling(rx-0.5) <= src_max_x-src_npts_bdr .and. &
 928           src_min_y+src_npts_bdr <= floor(ry+0.5) .and. ceiling(ry-0.5) <= src_max_y-src_npts_bdr) then
 929          is_point_in_tile = .true.
 930       else
 931          is_point_in_tile = .false.
 932       end if
 933
 934       return
 935
 936    end function is_point_in_tile
 937
 938 end module proc_point_module