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