geogrid/src/interp_module.F

   1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   2 ! MODULE INTERP_MODULE
   3 !
   4 ! This module provides routines for interpolation.
   5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   6 module interp_module
   7
   8    use bitarray_module
   9    use misc_definitions_module
  10    use module_debug
  11    use queue_module
  12
  13    contains
  14
  15    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  16    ! Name: interp_array_from_string
  17    !
  18    ! Purpose:
  19    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  20    function interp_array_from_string(interp_string)
  21
  22       implicit none
  23
  24       ! Arguments
  25       character (len=*), intent(in) :: interp_string
  26
  27       ! Local variables
  28       integer :: j, p1, p2, iend, num_methods
  29
  30       ! Return value
  31       integer, pointer, dimension(:) :: interp_array_from_string
  32
  33       ! Get an idea of how many interpolation methods are in the string
  34       !   so we can allocate an appropriately sized array
  35       num_methods = 1
  36       do j=1,len_trim(interp_string)
  37          if (interp_string(j:j) == '+') num_methods = num_methods + 1
  38       end do
  39
  40       allocate(interp_array_from_string(num_methods+1))
  41       interp_array_from_string = 0
  42
  43       iend = len_trim(interp_string)
  44
  45       p1 = 1
  46       p2 = index(interp_string(1:iend),'+')
  47       j = 1
  48       do while(p2 >= p1)
  49          if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
  50              len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
  51             interp_array_from_string(j) = N_NEIGHBOR
  52             j = j + 1
  53          else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
  54              len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
  55             interp_array_from_string(j) = AVERAGE4
  56             j = j + 1
  57          else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
  58              len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
  59             interp_array_from_string(j) = AVERAGE16
  60             j = j + 1
  61          else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
  62              len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
  63             interp_array_from_string(j) = W_AVERAGE4
  64             j = j + 1
  65          else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
  66              len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
  67             interp_array_from_string(j) = W_AVERAGE16
  68             j = j + 1
  69          else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
  70              len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
  71             interp_array_from_string(j) = FOUR_POINT
  72             j = j + 1
  73          else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
  74              len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
  75             interp_array_from_string(j) = SIXTEEN_POINT
  76             j = j + 1
  77          else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
  78              len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
  79             interp_array_from_string(j) = SEARCH
  80             j = j + 1
  81          else
  82             if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
  83                call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
  84          end if
  85          p1 = p2 + 1
  86          p2 = index(interp_string(p1:iend),'+') + p1 - 1
  87       end do
  88
  89       p2 = iend+1
  90       if (p1 < iend) then
  91          if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
  92              len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
  93             interp_array_from_string(j) = N_NEIGHBOR
  94             j = j + 1
  95          else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
  96              len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
  97             interp_array_from_string(j) = AVERAGE4
  98             j = j + 1
  99          else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
 100              len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
 101             interp_array_from_string(j) = AVERAGE16
 102             j = j + 1
 103          else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
 104              len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
 105             interp_array_from_string(j) = W_AVERAGE4
 106             j = j + 1
 107          else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
 108              len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
 109             interp_array_from_string(j) = W_AVERAGE16
 110             j = j + 1
 111          else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
 112              len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
 113             interp_array_from_string(j) = FOUR_POINT
 114             j = j + 1
 115          else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
 116              len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
 117             interp_array_from_string(j) = SIXTEEN_POINT
 118             j = j + 1
 119          else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
 120              len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
 121             interp_array_from_string(j) = SEARCH
 122             j = j + 1
 123          else
 124             if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
 125                call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
 126          end if
 127       end if
 128
 129       return
 130
 131    end function interp_array_from_string
 132
 133
 134    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 135    ! Name: interp_sequence
 136    !
 137    ! Purpose: Delegates the actual task of interpolation to specific
 138    !   interpolation routines defined in the module.
 139    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 140    recursive function interp_sequence(xx, yy, izz, array, start_x, end_x, &
 141               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 142
 143       implicit none
 144
 145       ! Arguments
 146       integer, intent(in) :: start_x, start_y, start_z
 147       integer, intent(in) :: end_x, end_y, end_z
 148       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 149                                                 !   interpolate within
 150       real, intent(in) :: xx , yy               ! The location to interpolate to
 151       real, intent(in) :: msgval
 152       real, intent(in), optional :: maskval
 153       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
 154       integer, intent(in) :: idx
 155       integer, dimension(:), intent(in) :: interp_list
 156       real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
 157       character (len=1), intent(in), optional :: mask_relational
 158
 159       ! Return value
 160       real :: interp_sequence
 161
 162       ! No more interpolation methods to try
 163       if (interp_list(idx) == 0) then
 164          interp_sequence = msgval
 165          return
 166       end if
 167
 168       if (interp_list(idx) == FOUR_POINT) then
 169          interp_sequence = four_pt(xx, yy, izz, array, start_x, end_x, &
 170                                    start_y, end_y, start_z, end_z, &
 171                                    msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 172       else if (interp_list(idx) == AVERAGE4) then
 173          interp_sequence = four_pt_average(xx, yy, izz, array, start_x, end_x, &
 174                                            start_y, end_y, start_z, end_z, &
 175                                            msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 176       else if (interp_list(idx) == W_AVERAGE4) then
 177          interp_sequence = wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
 178                                               start_y, end_y, start_z, end_z, &
 179                                               msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 180       else if (interp_list(idx) == N_NEIGHBOR) then
 181          interp_sequence = nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
 182                                             start_y, end_y, start_z, end_z, &
 183                                             msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 184       else if (interp_list(idx) == SIXTEEN_POINT) then
 185          interp_sequence = sixteen_pt(xx, yy, izz, array, start_x, end_x, &
 186                                       start_y, end_y, start_z, end_z, &
 187                                       msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 188       else if (interp_list(idx) == SEARCH) then
 189          interp_sequence = search_extrap(xx, yy, izz, array, start_x, end_x, &
 190                                          start_y, end_y, start_z, end_z, &
 191                                          msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 192       else if (interp_list(idx) == AVERAGE16) then
 193          interp_sequence = sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
 194                                               start_y, end_y, start_z, end_z, &
 195                                               msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 196       else if (interp_list(idx) == W_AVERAGE16) then
 197          interp_sequence = wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
 198                                                  start_y, end_y, start_z, end_z, &
 199                                                  msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
 200       end if
 201
 202    end function interp_sequence
 203
 204
 205    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 206    ! Name: nearest_neighbor
 207    !
 208    ! Purpose: Returns the point nearest to (xx,yy). If (xx,yy) is outside of the
 209    !   array, the point on the edge of the array nearest to (xx,yy) is returned.
 210    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 211    recursive function nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
 212               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 213
 214       implicit none
 215
 216       ! Arguments
 217       integer, intent(in) :: start_x, start_y, start_z
 218       integer, intent(in) :: end_x, end_y, end_z
 219       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 220                                                 !   interpolate within
 221       real, intent(in) :: xx , yy               ! The location to interpolate to
 222       real, intent(in) :: msgval
 223       real, intent(in), optional :: maskval
 224       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
 225       integer, dimension(:), intent(in) :: interp_list
 226       integer, intent(in) :: idx
 227       real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
 228       character (len=1), intent(in), optional :: mask_relational
 229
 230       ! Return value
 231       real :: nearest_neighbor
 232
 233       ! Local variables
 234       integer :: ix, iy
 235
 236       ix = nint(xx)
 237       iy = nint(yy)
 238
 239       ! The first thing to do is to ensure that the point (xx,yy) is within the array
 240       if (ix < start_x .or. ix > end_x) then
 241          nearest_neighbor = msgval
 242          return
 243       end if
 244
 245       if (iy < start_y .or. iy > end_y) then
 246          nearest_neighbor = msgval
 247          return
 248       end if
 249
 250       if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 251          if (mask_relational == '<' .and. mask_array(ix,iy) < maskval) then
 252             nearest_neighbor = msgval
 253          else if (mask_relational == '>' .and. mask_array(ix,iy) > maskval) then
 254             nearest_neighbor = msgval
 255          else if (mask_relational == ' ' .and. mask_array(ix,iy) == maskval) then
 256             nearest_neighbor = msgval
 257          else
 258             nearest_neighbor = array(ix,iy,izz)
 259          end if
 260       else if (present(mask_array) .and. present(maskval)) then
 261          if (maskval == mask_array(ix,iy)) then
 262             nearest_neighbor = msgval
 263          else
 264             nearest_neighbor = array(ix,iy,izz)
 265          end if
 266       else
 267          nearest_neighbor = array(ix,iy,izz)
 268       end if
 269
 270       ! If we have a missing value, try the next interpolation method in the sequence
 271       if (nearest_neighbor == msgval) then
 272          nearest_neighbor = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 273                              start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 274       end if
 275
 276    end function nearest_neighbor
 277
 278
 279    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 280    ! Name: search_extrap
 281    !
 282    ! Purpose: Returns the point nearest to (xx,yy) that has a non-missing value.
 283    !   If no valid value can be found in the array, msgval is returned.
 284    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 285    recursive function search_extrap(xx, yy, izz, array, start_x, end_x, &
 286               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 287
 288       implicit none
 289
 290       ! Arguments
 291       integer, intent(in) :: start_x, start_y, start_z
 292       integer, intent(in) :: end_x, end_y, end_z
 293       integer, intent(in) :: izz                ! The z-index of the 2d-array to search
 294       real, intent(in) :: xx , yy               ! The location of the search origin
 295       real, intent(in) :: msgval
 296       real, intent(in), optional :: maskval
 297       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
 298           intent(in) :: array
 299       integer, dimension(:), intent(in) :: interp_list
 300       integer, intent(in) :: idx
 301       real, dimension(start_x:end_x, start_y:end_y), &
 302           intent(in), optional :: mask_array
 303       character (len=1), intent(in), optional :: mask_relational
 304
 305       ! Return value
 306       real :: search_extrap
 307
 308       ! Local variables
 309       integer :: i, j
 310       real :: distance
 311       logical :: found_valid
 312       type (bitarray) :: b
 313       type (queue) :: q
 314       type (q_data) :: qdata
 315
 316       ! We only search if the starting point is within the array
 317       if (nint(xx) < start_x .or. nint(xx) > end_x .or. &
 318           nint(yy) < start_y .or. nint(yy) > end_y) then
 319          search_extrap = msgval
 320          return
 321       end if
 322
 323       call bitarray_create(b, (end_x-start_x+1), (end_y-start_y+1))
 324       call q_init(q)
 325
 326       found_valid = .false.
 327       qdata%x = nint(xx)
 328       qdata%y = nint(yy)
 329       call q_insert(q, qdata)
 330       call bitarray_set(b, qdata%x-start_x+1, qdata%y-start_y+1)
 331
 332       do while (q_isdata(q) .and. (.not. found_valid))
 333          qdata = q_remove(q)
 334          i = qdata%x
 335          j = qdata%y
 336
 337          if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 338             if (mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
 339                found_valid = .true.
 340             else if (mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
 341                found_valid = .true.
 342             else if (mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
 343                found_valid = .true.
 344             end if
 345          else if (present(mask_array) .and. present(maskval)) then
 346             if (array(i,j,izz) /= msgval .and. mask_array(i,j) /= maskval) found_valid = .true.
 347          else
 348             if (array(i,j,izz) /= msgval) found_valid = .true.
 349          end if
 350
 351          if (i-1 >= start_x) then
 352             if (.not. bitarray_test(b, (i-1)-start_x+1, j-start_y+1)) then
 353                qdata%x = i-1
 354                qdata%y = j
 355                call q_insert(q, qdata)
 356                call bitarray_set(b, (i-1)-start_x+1, j-start_y+1)
 357             end if
 358          end if
 359          if (i+1 <= end_x) then
 360             if (.not. bitarray_test(b, (i+1)-start_x+1, j-start_y+1)) then
 361                qdata%x = i+1
 362                qdata%y = j
 363                call q_insert(q, qdata)
 364                call bitarray_set(b, (i+1)-start_x+1, j-start_y+1)
 365             end if
 366          end if
 367          if (j-1 >= start_y) then
 368             if (.not. bitarray_test(b, i-start_x+1, (j-1)-start_y+1)) then
 369                qdata%x = i
 370                qdata%y = j-1
 371                call q_insert(q, qdata)
 372                call bitarray_set(b, i-start_x+1, (j-1)-start_y+1)
 373             end if
 374          end if
 375          if (j+1 <= end_y) then
 376             if (.not. bitarray_test(b, i-start_x+1, (j+1)-start_y+1)) then
 377                qdata%x = i
 378                qdata%y = j+1
 379                call q_insert(q, qdata)
 380                call bitarray_set(b, i-start_x+1, (j+1)-start_y+1)
 381             end if
 382          end if
 383       end do
 384
 385       if (found_valid) then
 386          distance = (real(i)-xx)*(real(i)-xx)+(real(j)-yy)*(real(j)-yy)
 387          search_extrap = array(i,j,izz)
 388          do while (q_isdata(q))
 389             qdata = q_remove(q)
 390             if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 391                if (mask_relational == '<' .and. mask_array(qdata%x,qdata%y) >= maskval &
 392                    .and. array(qdata%x,qdata%y,izz) /= msgval) then
 393                   if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
 394                      distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
 395                      search_extrap = array(qdata%x, qdata%y, izz)
 396                   end if
 397                else if (mask_relational == '>' .and. mask_array(qdata%x,qdata%y) <= maskval &
 398                         .and. array(qdata%x,qdata%y,izz) /= msgval) then
 399                   if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
 400                      distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
 401                      search_extrap = array(qdata%x, qdata%y, izz)
 402                   end if
 403                else if (mask_relational == ' ' .and. mask_array(qdata%x,qdata%y) /= maskval &
 404                         .and. array(qdata%x,qdata%y,izz) /= msgval) then
 405                   if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
 406                      distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
 407                      search_extrap = array(qdata%x, qdata%y, izz)
 408                   end if
 409                end if
 410
 411             else if (present(mask_array) .and. present(maskval)) then
 412                if (array(qdata%x,qdata%y,izz) /= msgval .and. mask_array(qdata%x,qdata%y) /= maskval) then
 413                   if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
 414                      distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
 415                      search_extrap = array(qdata%x, qdata%y, izz)
 416                   end if
 417                end if
 418
 419             else
 420                if (array(qdata%x,qdata%y,izz) /= msgval) then
 421                   if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
 422                      distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
 423                      search_extrap = array(qdata%x, qdata%y, izz)
 424                   end if
 425                end if
 426             end if
 427          end do
 428       else
 429          search_extrap = msgval
 430       end if
 431
 432       call q_destroy(q)
 433       call bitarray_destroy(b)
 434
 435    end function search_extrap
 436
 437
 438    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 439    ! Name: four_pt_average
 440    !
 441    ! Purpose: Average of four surrounding grid point values
 442    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 443    recursive function four_pt_average(xx, yy, izz, array, start_x, end_x, &
 444               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 445
 446       implicit none
 447
 448       ! Arguments
 449       integer, intent(in) :: start_x, start_y, start_z
 450       integer, intent(in) :: end_x, end_y, end_z
 451       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 452                                                 !   interpolate within
 453       real, intent(in) :: xx, yy                ! The location to interpolate to
 454       real, intent(in) :: msgval
 455       real, intent(in), optional :: maskval
 456       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
 457           intent(in) :: array
 458       integer, dimension(:), intent(in) :: interp_list
 459       integer, intent(in) :: idx
 460       real, dimension(start_x:end_x, start_y:end_y), &
 461           intent(in), optional :: mask_array
 462       character (len=1), intent(in), optional :: mask_relational
 463
 464       ! Return value
 465       real :: four_pt_average
 466
 467       ! Local variables
 468       integer :: ifx, ify, icx, icy
 469       real :: fxfy, fxcy, cxfy, cxcy
 470
 471       fxfy = 1.0
 472       fxcy = 1.0
 473       cxfy = 1.0
 474       cxcy = 1.0
 475
 476       ifx = floor(xx)
 477       icx = ceiling(xx)
 478       ify = floor(yy)
 479       icy = ceiling(yy)
 480
 481       ! First, make sure that the point is contained in the source array
 482       if (ifx < start_x .or. icx > end_x .or. &
 483           ify < start_y .or. icy > end_y) then
 484
 485          ! But if the point is at most half a grid point out, we can
 486          !   still proceed with modified ifx, icx, ify, and icy.
 487          if (xx > real(start_x)-0.5 .and. ifx < start_x) then
 488             ifx = start_x
 489             icx = start_x
 490          else if (xx < real(end_x)+0.5 .and. icx > end_x) then
 491             ifx = end_x
 492             icx = end_x
 493          end if
 494
 495          if (yy > real(start_y)-0.5 .and. ify < start_y) then
 496             ify = start_y
 497             icy = start_y
 498          else if (yy < real(end_y)+0.5 .and. icy > end_y) then
 499             ify = end_y
 500             icy = end_y
 501          end if
 502
 503          if (ifx < start_x .or. icx > end_x .or. &
 504              ify < start_y .or. icy > end_y) then
 505             four_pt_average = msgval
 506             return
 507          end if
 508       end if
 509
 510       if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 511          ! we determine which maskval is useable by... if the symbol > is found, then only
 512          !    values less than 'maskval' can be used and if the symbol < is found,
 513          !    then only the values greater than 'maskval' can be used.
 514          if (mask_relational == '>') then
 515             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0
 516             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0
 517             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0
 518             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0
 519          else if (mask_relational == '<') then
 520             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0
 521             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0
 522             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0
 523             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
 524          else  !equal
 525             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
 526             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
 527             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
 528             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
 529          end if
 530       else if (present(mask_array) .and. present(maskval)) then
 531          if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
 532          if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
 533          if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
 534          if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
 535       else
 536          if (array(ifx, ify, izz) == msgval) fxfy = 0.0
 537          if (array(ifx, icy, izz) == msgval) fxcy = 0.0
 538          if (array(icx, ify, izz) == msgval) cxfy = 0.0
 539          if (array(icx, icy, izz) == msgval) cxcy = 0.0
 540       end if
 541
 542       ! If all four points are missing, try the next interpolation method in the sequence
 543       if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
 544          four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 545                              start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 546       else
 547          four_pt_average = (fxfy * array(ifx, ify, izz) + &
 548                             fxcy * array(ifx, icy, izz) + &
 549                             cxfy * array(icx, ify, izz) + &
 550                             cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
 551       end if
 552
 553    end function four_pt_average
 554
 555
 556    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 557    ! Name: wt_four_pt_average
 558    !
 559    ! Purpose: Weighted average of four surrounding grid point values
 560    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 561    recursive function wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
 562               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 563
 564       implicit none
 565
 566       ! Arguments
 567       integer, intent(in) :: start_x, start_y, start_z
 568       integer, intent(in) :: end_x, end_y, end_z
 569       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 570                                                 !   interpolate within
 571       real, intent(in) :: xx, yy                ! The location to interpolate to
 572       real, intent(in) :: msgval
 573       real, intent(in), optional :: maskval
 574       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
 575           intent(in) :: array
 576       integer, dimension(:), intent(in) :: interp_list
 577       integer, intent(in) :: idx
 578       real, dimension(start_x:end_x, start_y:end_y), &
 579           intent(in), optional :: mask_array
 580       character (len=1), intent(in), optional :: mask_relational
 581
 582       ! Return value
 583       real :: wt_four_pt_average
 584
 585       ! Local variables
 586       integer :: ifx, ify, icx, icy
 587       real :: fxfy, fxcy, cxfy, cxcy
 588
 589       ifx = floor(xx)
 590       icx = ceiling(xx)
 591       ify = floor(yy)
 592       icy = ceiling(yy)
 593
 594       fxfy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(ify))**2))
 595       fxcy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(icy))**2))
 596       cxfy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(ify))**2))
 597       cxcy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(icy))**2))
 598
 599       ! First, make sure that the point is contained in the source array
 600       if (ifx < start_x .or. icx > end_x .or. &
 601           ify < start_y .or. icy > end_y) then
 602
 603          ! But if the point is at most half a grid point out, we can
 604          !   still proceed with modified ifx, icx, ify, and icy.
 605          if (xx > real(start_x)-0.5 .and. ifx < start_x) then
 606             ifx = start_x
 607             icx = start_x
 608          else if (xx < real(end_x)+0.5 .and. icx > end_x) then
 609             ifx = end_x
 610             icx = end_x
 611          end if
 612
 613          if (yy > real(start_y)-0.5 .and. ifx < start_y) then
 614             ify = start_y
 615             icy = start_y
 616          else if (yy < real(end_y)+0.5 .and. icy > end_y) then
 617             ify = end_y
 618             icy = end_y
 619          end if
 620
 621          if (ifx < start_x .or. icx > end_x .or. &
 622              ify < start_y .or. icy > end_y) then
 623             wt_four_pt_average = msgval
 624             return
 625          end if
 626       end if
 627
 628       if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 629          ! we determine which maskval is useable by... if the symbol > is found, then only
 630          !    values less than 'maskval' can be used and if the symbol < is found,
 631          !    then only the values greater than 'maskval' can be used.
 632          if (mask_relational == '>') then
 633             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0
 634             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0
 635             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0
 636             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0
 637          else if (mask_relational == '<') then
 638             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0
 639             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0
 640             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0
 641             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
 642          else  !equal
 643             if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
 644             if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
 645             if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
 646             if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
 647          end if
 648       else if (present(mask_array) .and. present(maskval)) then
 649          if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
 650          if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
 651          if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
 652          if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
 653       else
 654          if (array(ifx, ify, izz) == msgval) fxfy = 0.0
 655          if (array(ifx, icy, izz) == msgval) fxcy = 0.0
 656          if (array(icx, ify, izz) == msgval) cxfy = 0.0
 657          if (array(icx, icy, izz) == msgval) cxcy = 0.0
 658       end if
 659
 660       ! If all four points are missing, try the next interpolation method in the sequence
 661       if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
 662          wt_four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 663                                 start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 664       else
 665          wt_four_pt_average = (fxfy * array(ifx, ify, izz) + &
 666                                fxcy * array(ifx, icy, izz) + &
 667                                cxfy * array(icx, ify, izz) + &
 668                                cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
 669       end if
 670
 671    end function wt_four_pt_average
 672
 673
 674    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 675    ! Name: sixteen_pt_average
 676    !
 677    ! Purpose: Average of sixteen surrounding grid point values
 678    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 679    recursive function sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
 680               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 681
 682       implicit none
 683
 684       ! Arguments
 685       integer, intent(in) :: start_x, start_y, start_z
 686       integer, intent(in) :: end_x, end_y, end_z
 687       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 688                                                 !   interpolate within
 689       real, intent(in) :: xx , yy               ! The location to interpolate to
 690       real, intent(in) :: msgval
 691       real, intent(in), optional :: maskval
 692       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
 693           intent(in) :: array
 694       integer, dimension(:), intent(in) :: interp_list
 695       integer, intent(in) :: idx
 696       real, dimension(start_x:end_x, start_y:end_y), &
 697           intent(in), optional :: mask_array
 698       character (len=1), intent(in), optional :: mask_relational
 699
 700       ! Return value
 701       real :: sixteen_pt_average
 702
 703       ! Local variables
 704       integer :: i, j, ifx, ify
 705       real :: sum, sum_weight
 706       real, dimension(4,4) :: weights
 707
 708       ifx = floor(xx)
 709       ify = floor(yy)
 710
 711       ! First see whether the point is far enough within the array to
 712       !   allow for a sixteen point average.
 713       if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
 714           ify < start_y+1 .or. ify > end_y-2) then
 715          sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 716                                 start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 717          return
 718       end if
 719
 720       sum_weight = 0.0
 721       do i=1,4
 722          do j=1,4
 723             if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 724                if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
 725                    weights(i,j) = 0.0
 726                else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
 727                    weights(i,j) = 0.0
 728                else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
 729                    weights(i,j) = 0.0
 730                else
 731                    weights(i,j) = 1.0
 732                end if
 733                if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
 734             else if (present(mask_array) .and. present(maskval)) then
 735                if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
 736                   weights(i,j) = 0.0
 737                else
 738                   weights(i,j) = 1.0
 739                end if
 740             else
 741                if (array(ifx+3-i, ify+3-j, izz) == msgval) then
 742                   weights(i,j) = 0.0
 743                else
 744                   weights(i,j) = 1.0
 745                end if
 746             end if
 747
 748             sum_weight = sum_weight + weights(i,j)
 749
 750          end do
 751       end do
 752
 753       ! If all points are missing, try the next interpolation method in the sequence
 754       if (sum_weight == 0.0) then
 755          sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 756                                 start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 757       else
 758          sum = 0.0
 759          do i=1,4
 760             do j=1,4
 761                sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
 762             end do
 763          end do
 764          sixteen_pt_average = sum / sum_weight
 765       end if
 766
 767    end function sixteen_pt_average
 768
 769
 770    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 771    ! Name: wt_sixteen_pt_average
 772    !
 773    ! Purpose: Weighted average of sixteen surrounding grid point values
 774    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 775    recursive function wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
 776               start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 777
 778       implicit none
 779
 780       ! Arguments
 781       integer, intent(in) :: start_x, start_y, start_z
 782       integer, intent(in) :: end_x, end_y, end_z
 783       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 784                                                 !   interpolate within
 785       real, intent(in) :: xx , yy               ! The location to interpolate to
 786       real, intent(in) :: msgval
 787       real, intent(in), optional :: maskval
 788       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
 789           intent(in) :: array
 790       integer, dimension(:), intent(in) :: interp_list
 791       integer, intent(in) :: idx
 792       real, dimension(start_x:end_x, start_y:end_y), &
 793           intent(in), optional :: mask_array
 794       character (len=1), intent(in), optional :: mask_relational
 795
 796       ! Return value
 797       real :: wt_sixteen_pt_average
 798
 799       ! Local variables
 800       integer :: i, j, ifx, ify
 801       real :: sum, sum_weight
 802       real, dimension(4,4) :: weights
 803
 804       ifx = floor(xx)
 805       ify = floor(yy)
 806
 807       ! First see whether the point is far enough within the array to
 808       !   allow for a sixteen point average.
 809       if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
 810           ify < start_y+1 .or. ify > end_y-2) then
 811          wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 812                                    start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 813          return
 814       end if
 815
 816       sum_weight = 0.0
 817       do i=1,4
 818          do j=1,4
 819             if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 820                if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
 821                    weights(i,j) = 0.0
 822                else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
 823                    weights(i,j) = 0.0
 824                else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
 825                    weights(i,j) = 0.0
 826                else
 827                    weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
 828                end if
 829                if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
 830             else if (present(mask_array) .and. present(maskval)) then
 831                if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
 832                   weights(i,j) = 0.0
 833                else
 834                   weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
 835                end if
 836             else
 837                if (array(ifx+3-i, ify+3-j, izz) == msgval) then
 838                   weights(i,j) = 0.0
 839                else
 840                   weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
 841                end if
 842             end if
 843
 844             sum_weight = sum_weight + weights(i,j)
 845
 846          end do
 847       end do
 848
 849       ! If all points are missing, try the next interpolation method in the sequence
 850       if (sum_weight == 0.0) then
 851          wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 852                                    start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 853       else
 854          sum = 0.0
 855          do i=1,4
 856             do j=1,4
 857                sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
 858             end do
 859          end do
 860          wt_sixteen_pt_average = sum / sum_weight
 861       end if
 862
 863    end function wt_sixteen_pt_average
 864
 865
 866    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 867    ! Name: four_pt
 868    !
 869    ! Purpose: Bilinear interpolation among four grid values
 870    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 871    recursive function four_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 872                        start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 873
 874       implicit none
 875
 876       ! Arguments
 877       integer, intent(in) :: start_x, start_y, start_z
 878       integer, intent(in) :: end_x, end_y, end_z
 879       integer, intent(in) :: izz                ! The z-index of the 2d-array to
 880                                                 !   interpolate within
 881       real, intent(in) :: xx , yy               ! The location to interpolate to
 882       real, intent(in) :: msgval
 883       real, intent(in), optional :: maskval
 884       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
 885       integer, dimension(:), intent(in) :: interp_list
 886       integer, intent(in) :: idx
 887       real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
 888       character (len=1), intent(in), optional :: mask_relational
 889
 890       ! Return value
 891       real :: four_pt
 892
 893       ! Local variables
 894       integer :: min_x, min_y, max_x, max_y
 895
 896       min_x = floor(xx)
 897       min_y = floor(yy)
 898       max_x = ceiling(xx)
 899       max_y = ceiling(yy)
 900
 901       if (min_x < start_x .or. max_x > end_x) then
 902          four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 903                                    msgval, interp_list, idx, mask_relational, maskval, mask_array)
 904          return
 905       end if
 906
 907       if (min_y < start_y .or. max_y > end_y) then
 908          four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 909                                    msgval, interp_list, idx, mask_relational, maskval, mask_array)
 910          return
 911       end if
 912
 913       ! If we have a missing value, try the next interpolation method in the sequence
 914       if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
 915          ! we determine which maskval is useable by... if the symbol > is found, then only
 916          !    values less than 'maskval' can be used and if the symbol < is found,
 917          !    then only the values greater than 'maskval' can be used.
 918          if (mask_relational == '>' ) then
 919             if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) > maskval .or. &
 920                 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) > maskval .or. &
 921                 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) > maskval .or. &
 922                 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) > maskval) then
 923                four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 924                                      msgval, interp_list, idx, mask_relational, maskval, mask_array)
 925                return
 926             end if
 927          else if (mask_relational == '<' ) then
 928             if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) < maskval .or. &
 929                 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) < maskval .or. &
 930                 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) < maskval .or. &
 931                 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) < maskval) then
 932                four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 933                                          msgval, interp_list, idx, mask_relational, maskval, mask_array)
 934                return
 935              end if
 936          else if (mask_relational == ' ' ) then
 937             if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
 938                 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
 939                 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
 940                 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then
 941                four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 942                                          msgval, interp_list, idx, mask_relational, maskval, mask_array)
 943                return
 944             end if
 945          end if
 946       else if (present(mask_array) .and. present(maskval)) then
 947          if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
 948              array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
 949              array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
 950              array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then
 951             four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 952                                       msgval, interp_list, idx, mask_relational, maskval, mask_array)
 953          end if
 954       else
 955          if (array(min_x,min_y,izz) == msgval .or. &
 956              array(max_x,min_y,izz) == msgval .or. &
 957              array(min_x,max_y,izz) == msgval .or. &
 958              array(max_x,max_y,izz) == msgval ) then
 959             four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
 960                                 start_z, end_z, msgval, interp_list, idx)
 961             return
 962          end if
 963       end if
 964
 965       if (min_x == max_x) then
 966          if (min_y == max_y) then
 967             four_pt = array(min_x,min_y,izz)
 968          else
 969             four_pt = array(min_x,min_y,izz)*(real(max_y)-yy) + &
 970                       array(min_x,max_y,izz)*(yy-real(min_y))
 971          end if
 972       else if (min_y == max_y) then
 973          if (min_x == max_x) then
 974             four_pt = array(min_x,min_y,izz)
 975          else
 976             four_pt = array(min_x,min_y,izz)*(real(max_x)-xx) + &
 977                       array(max_x,min_y,izz)*(xx-real(min_x))
 978          end if
 979       else
 980          four_pt = (yy - min_y) * (array(min_x,max_y,izz)*(real(max_x)-xx) + &
 981                                    array(max_x,max_y,izz)*(xx-real(min_x))) + &
 982                    (max_y - yy) * (array(min_x,min_y,izz)*(real(max_x)-xx) + &
 983                                    array(max_x,min_y,izz)*(xx-real(min_x)));
 984       end if
 985
 986    end function four_pt
 987
 988
 989    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 990    ! Name: sixteen_pt
 991    !
 992    ! Purpose: Overlapping parabolic interpolation among sixteen grid values
 993    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 994    recursive function sixteen_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
 995                                  msgval, interp_list, idx, mask_relational, maskval, mask_array)
 996
 997       implicit none
 998
 999       ! Arguments
1000       integer, intent(in) :: izz    ! z-index of 2d-array to interpolate within
1001       integer, intent(in) :: start_x, start_y, start_z
1002       integer, intent(in) :: end_x, end_y, end_z
1003       real, intent(in) :: xx , yy              ! The location to interpolate to
1004       real, intent(in) :: msgval
1005       real, intent(in), optional :: maskval
1006       real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
1007           intent(in) :: array
1008       integer, dimension(:), intent(in) :: interp_list
1009       integer, intent(in) :: idx
1010       real, dimension(start_x:end_x, start_y:end_y), &
1011           intent(in), optional :: mask_array
1012       character (len=1), intent(in), optional :: mask_relational
1013
1014       ! Return value
1015       real :: sixteen_pt
1016
1017       ! Local variables
1018       integer :: n , i , j , k , kk , l , ll
1019       real :: x , y , a , b , c , d , e , f , g , h
1020       real, dimension(4,4) :: stl
1021       logical :: is_masked
1022
1023       is_masked = .false.
1024
1025       if (int(xx) < start_x .or. int(xx) > end_x .or. &
1026           int(yy) < start_y .or. int(yy) > end_y) then
1027          sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1028                                       start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
1029          return
1030       end if
1031
1032       sixteen_pt = 0.0
1033       n = 0
1034       i = int(xx + 0.00001)
1035       j = int(yy + 0.00001)
1036       x = xx - i
1037       y = yy - j
1038
1039       if ( ( abs(x) > 0.0001 ) .or. ( abs(y) > 0.0001 ) ) then
1040
1041          loop_1 : do k = 1,4
1042             kk = i + k - 2
1043             if ( kk < start_x) then
1044                kk = start_x
1045             else if ( kk > end_x) then
1046                kk = end_x
1047             end if
1048             loop_2 : do l = 1,4
1049                stl(k,l) = 0.
1050                ll = j + l - 2
1051                if ( ll < start_y ) then
1052                   ll = start_y
1053                else if ( ll > end_y) then
1054                   ll = end_y
1055                end if
1056                stl(k,l) = array(kk,ll,izz)
1057                n = n + 1
1058                if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1059                   if (mask_relational == '>' .and. mask_array(kk,ll) > maskval) then
1060                      is_masked = .true.
1061                   else if (mask_relational == '<' .and. mask_array(kk,ll) < maskval) then
1062                      is_masked = .true.
1063                   else if (mask_relational == ' ' .and. mask_array(kk,ll) == maskval) then
1064                      is_masked = .true.
1065                   end if
1066                else if (present(mask_array) .and. present(maskval)) then
1067                   if (mask_array(kk,ll) == maskval) is_masked = .true.
1068                end if
1069                if ( stl(k,l) == 0. .and. msgval /= 0.) then
1070                   stl(k,l) = 1.E-20
1071                end if
1072             end do loop_2
1073          end do loop_1
1074
1075          ! If we have a missing value, try the next interpolation method in the sequence
1076          if (present(mask_array) .and. present(maskval)) then
1077             do k=1,4
1078                do l=1,4
1079                   if (stl(k,l) == msgval .or. is_masked) then
1080                      sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1081                                                   start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
1082                      return
1083                   end if
1084                end do
1085             end do
1086          else
1087             do k=1,4
1088                do l=1,4
1089                   if (stl(k,l) == msgval) then
1090                      sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1091                                                   start_z, end_z, msgval, interp_list, idx)
1092                      return
1093                   end if
1094                end do
1095             end do
1096          end if
1097
1098          a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
1099          b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
1100          c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
1101          d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
1102          sixteen_pt = oned(y,a,b,c,d)
1103
1104          if (n /= 16) then
1105             e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
1106             f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
1107             g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
1108             h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
1109             sixteen_pt = (sixteen_pt+oned(x,e,f,g,h)) * 0.5
1110          end if
1111
1112          if (sixteen_pt == 1.E-20) sixteen_pt = 0.
1113
1114       else
1115          if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1116             if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1117                mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
1118                sixteen_pt = array(i,j,izz)
1119             else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1120                mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
1121                sixteen_pt = array(i,j,izz)
1122             else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1123                mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
1124                sixteen_pt = array(i,j,izz)
1125             else
1126                sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1127                                             msgval, interp_list, idx, mask_relational, maskval, mask_array)
1128             end if
1129          else if (present(mask_array) .and. present(maskval)) then
1130             if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1131                 mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
1132                sixteen_pt = array(i,j,izz)
1133             else
1134                sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1135                                             msgval, interp_list, idx, mask_relational, maskval, mask_array)
1136             end if
1137          else
1138             if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. array(i,j,izz) /= msgval) then
1139                sixteen_pt = array(i,j,izz)
1140             else
1141                sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1142                                             msgval, interp_list, idx, mask_relational, maskval, mask_array)
1143             end if
1144          end if
1145       end if
1146
1147    end function sixteen_pt
1148
1149
1150    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1151    ! Name: oned
1152    !
1153    ! Purpose: 1-dimensional overlapping parabolic interpolation
1154    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1155    function oned(x,a,b,c,d)
1156
1157       implicit none
1158
1159       ! Arguments
1160       real, intent(in) :: x,a,b,c,d
1161
1162       ! Return value
1163       real :: oned
1164
1165       oned = 0.
1166
1167       if ( x == 0. ) then
1168          oned = b
1169       else if ( x == 1. ) then
1170          oned = c
1171       end if
1172
1173       if (b*c /= 0.) then
1174          if ( a*d == 0. ) then
1175             if ( ( a == 0 ) .and. ( d == 0 ) ) then
1176                oned = b*(1.0-x)+c*x
1177             else if ( a /= 0. ) then
1178                oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
1179             else if ( d /= 0. ) then
1180                oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
1181             end if
1182          else
1183             oned = (1.0-x)*(b+x*(0.5*(c-a)+x*(0.5*(c+a)-b)))+x*(c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)))
1184          end if
1185       end if
1186
1187    end function oned
1188
1189 end module interp_module