/WPS/geogrid/src/interp_module.F

Large files files are truncated, but you can click here to view the full file

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! MODULE INTERP_MODULE
!
! This module provides routines for interpolation.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module interp_module

   use bitarray_module
   use misc_definitions_module
   use module_debug
   use queue_module
 
   contains

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: interp_array_from_string
   !
   ! Purpose:
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   function interp_array_from_string(interp_string)

      implicit none

      ! Arguments
      character (len=*), intent(in) :: interp_string

      ! Local variables
      integer :: j, p1, p2, iend, num_methods 

      ! Return value
      integer, pointer, dimension(:) :: interp_array_from_string

      ! Get an idea of how many interpolation methods are in the string
      !   so we can allocate an appropriately sized array
      num_methods = 1
      do j=1,len_trim(interp_string)
         if (interp_string(j:j) == '+') num_methods = num_methods + 1
      end do

      allocate(interp_array_from_string(num_methods+1))
      interp_array_from_string = 0

      iend = len_trim(interp_string)

      p1 = 1
      p2 = index(interp_string(1:iend),'+')
      j = 1
      do while(p2 >= p1)
         if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
            interp_array_from_string(j) = N_NEIGHBOR
            j = j + 1
         else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
            interp_array_from_string(j) = AVERAGE4
            j = j + 1
         else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
            interp_array_from_string(j) = AVERAGE16
            j = j + 1
         else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
            interp_array_from_string(j) = W_AVERAGE4
            j = j + 1
         else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
            interp_array_from_string(j) = W_AVERAGE16
            j = j + 1
         else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
            interp_array_from_string(j) = FOUR_POINT
            j = j + 1
         else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
            interp_array_from_string(j) = SIXTEEN_POINT
            j = j + 1
         else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
            interp_array_from_string(j) = SEARCH
            j = j + 1
         else
            if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
               call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
         end if
         p1 = p2 + 1
         p2 = index(interp_string(p1:iend),'+') + p1 - 1
      end do

      p2 = iend+1
      if (p1 < iend) then
         if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
            interp_array_from_string(j) = N_NEIGHBOR
            j = j + 1
         else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
            interp_array_from_string(j) = AVERAGE4
            j = j + 1
         else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
            interp_array_from_string(j) = AVERAGE16
            j = j + 1
         else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
            interp_array_from_string(j) = W_AVERAGE4
            j = j + 1
         else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
            interp_array_from_string(j) = W_AVERAGE16
            j = j + 1
         else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
            interp_array_from_string(j) = FOUR_POINT
            j = j + 1
         else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
            interp_array_from_string(j) = SIXTEEN_POINT
            j = j + 1
         else if (index(interp_string(p1:p2-1),'search') /= 0 .and. &
             len_trim(interp_string(p1:p2-1)) == len_trim('search')) then
            interp_array_from_string(j) = SEARCH
            j = j + 1
         else
            if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
               call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
         end if
      end if

      return

   end function interp_array_from_string
 

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: interp_sequence
   !
   ! Purpose: Delegates the actual task of interpolation to specific 
   !   interpolation routines defined in the module.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function interp_sequence(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to 
                                                !   interpolate within
      real, intent(in) :: xx , yy               ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
      integer, intent(in) :: idx
      integer, dimension(:), intent(in) :: interp_list
      real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational
  
      ! Return value
      real :: interp_sequence       
  
      ! No more interpolation methods to try
      if (interp_list(idx) == 0) then
         interp_sequence = msgval
         return
      end if
  
      if (interp_list(idx) == FOUR_POINT) then
         interp_sequence = four_pt(xx, yy, izz, array, start_x, end_x, &
                                   start_y, end_y, start_z, end_z, &
                                   msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == AVERAGE4) then
         interp_sequence = four_pt_average(xx, yy, izz, array, start_x, end_x, &
                                           start_y, end_y, start_z, end_z, &
                                           msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == W_AVERAGE4) then
         interp_sequence = wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
                                              start_y, end_y, start_z, end_z, &
                                              msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == N_NEIGHBOR) then
         interp_sequence = nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
                                            start_y, end_y, start_z, end_z, &
                                            msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == SIXTEEN_POINT) then
         interp_sequence = sixteen_pt(xx, yy, izz, array, start_x, end_x, &
                                      start_y, end_y, start_z, end_z, &
                                      msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == SEARCH) then
         interp_sequence = search_extrap(xx, yy, izz, array, start_x, end_x, &
                                         start_y, end_y, start_z, end_z, &
                                         msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == AVERAGE16) then
         interp_sequence = sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
                                              start_y, end_y, start_z, end_z, &
                                              msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      else if (interp_list(idx) == W_AVERAGE16) then
         interp_sequence = wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
                                                 start_y, end_y, start_z, end_z, &
                                                 msgval, interp_list, idx+1, mask_relational, maskval, mask_array)
      end if
 
   end function interp_sequence
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: nearest_neighbor
   !
   ! Purpose: Returns the point nearest to (xx,yy). If (xx,yy) is outside of the
   !   array, the point on the edge of the array nearest to (xx,yy) is returned.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
    
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to 
                                                !   interpolate within
      real, intent(in) :: xx , yy               ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: nearest_neighbor
  
      ! Local variables
      integer :: ix, iy

      ix = nint(xx)
      iy = nint(yy)
  
      ! The first thing to do is to ensure that the point (xx,yy) is within the array
      if (ix < start_x .or. ix > end_x) then
         nearest_neighbor = msgval
         return
      end if
  
      if (iy < start_y .or. iy > end_y) then
         nearest_neighbor = msgval
         return
      end if

      if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
         if (mask_relational == '<' .and. mask_array(ix,iy) < maskval) then 
            nearest_neighbor = msgval
         else if (mask_relational == '>' .and. mask_array(ix,iy) > maskval) then 
            nearest_neighbor = msgval
         else if (mask_relational == ' ' .and. mask_array(ix,iy) == maskval) then
            nearest_neighbor = msgval
         else
            nearest_neighbor = array(ix,iy,izz)
         end if        
      else if (present(mask_array) .and. present(maskval)) then 
         if (maskval == mask_array(ix,iy)) then
            nearest_neighbor = msgval
         else
            nearest_neighbor = array(ix,iy,izz)
         end if        
      else
         nearest_neighbor = array(ix,iy,izz)
      end if
  
      ! If we have a missing value, try the next interpolation method in the sequence
      if (nearest_neighbor == msgval) then
         nearest_neighbor = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                             start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
      end if
  
   end function nearest_neighbor
 

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: search_extrap
   !
   ! Purpose: Returns the point nearest to (xx,yy) that has a non-missing value.
   !   If no valid value can be found in the array, msgval is returned.
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function search_extrap(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to search
      real, intent(in) :: xx , yy               ! The location of the search origin 
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: search_extrap
  
      ! Local variables
      integer :: i, j
      real :: distance
      logical :: found_valid
      type (bitarray) :: b
      type (queue) :: q
      type (q_data) :: qdata

      ! We only search if the starting point is within the array
      if (nint(xx) < start_x .or. nint(xx) > end_x .or. &
          nint(yy) < start_y .or. nint(yy) > end_y) then
         search_extrap = msgval
         return 
      end if
      
      call bitarray_create(b, (end_x-start_x+1), (end_y-start_y+1))
      call q_init(q)
  
      found_valid = .false.
      qdata%x = nint(xx)
      qdata%y = nint(yy)
      call q_insert(q, qdata)
      call bitarray_set(b, qdata%x-start_x+1, qdata%y-start_y+1)
  
      do while (q_isdata(q) .and. (.not. found_valid))
         qdata = q_remove(q) 
         i = qdata%x
         j = qdata%y

         if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
            if (mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
               found_valid = .true.
            else if (mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
               found_valid = .true.
            else if (mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
               found_valid = .true.
            end if
         else if (present(mask_array) .and. present(maskval)) then
            if (array(i,j,izz) /= msgval .and. mask_array(i,j) /= maskval) found_valid = .true.                  
         else
            if (array(i,j,izz) /= msgval) found_valid = .true. 
         end if

         if (i-1 >= start_x) then
            if (.not. bitarray_test(b, (i-1)-start_x+1, j-start_y+1)) then
               qdata%x = i-1
               qdata%y = j
               call q_insert(q, qdata)
               call bitarray_set(b, (i-1)-start_x+1, j-start_y+1)
            end if
         end if  
         if (i+1 <= end_x) then
            if (.not. bitarray_test(b, (i+1)-start_x+1, j-start_y+1)) then
               qdata%x = i+1
               qdata%y = j
               call q_insert(q, qdata)
               call bitarray_set(b, (i+1)-start_x+1, j-start_y+1)
            end if
         end if  
         if (j-1 >= start_y) then
            if (.not. bitarray_test(b, i-start_x+1, (j-1)-start_y+1)) then
               qdata%x = i
               qdata%y = j-1
               call q_insert(q, qdata)
               call bitarray_set(b, i-start_x+1, (j-1)-start_y+1)
            end if
         end if  
         if (j+1 <= end_y) then
            if (.not. bitarray_test(b, i-start_x+1, (j+1)-start_y+1)) then
               qdata%x = i
               qdata%y = j+1
               call q_insert(q, qdata)
               call bitarray_set(b, i-start_x+1, (j+1)-start_y+1)
            end if
         end if  
      end do
  
      if (found_valid) then
         distance = (real(i)-xx)*(real(i)-xx)+(real(j)-yy)*(real(j)-yy) 
         search_extrap = array(i,j,izz)
         do while (q_isdata(q))
            qdata = q_remove(q)
            if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
               if (mask_relational == '<' .and. mask_array(qdata%x,qdata%y) >= maskval &
                   .and. array(qdata%x,qdata%y,izz) /= msgval) then
                  if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
                     distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) 
                     search_extrap = array(qdata%x, qdata%y, izz)
                  end if
               else if (mask_relational == '>' .and. mask_array(qdata%x,qdata%y) <= maskval &
                        .and. array(qdata%x,qdata%y,izz) /= msgval) then
                  if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
                     distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) 
                     search_extrap = array(qdata%x, qdata%y, izz)
                  end if
               else if (mask_relational == ' ' .and. mask_array(qdata%x,qdata%y) /= maskval &
                        .and. array(qdata%x,qdata%y,izz) /= msgval) then
                  if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
                     distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) 
                     search_extrap = array(qdata%x, qdata%y, izz)
                  end if
               end if
                
            else if (present(mask_array) .and. present(maskval)) then
               if (array(qdata%x,qdata%y,izz) /= msgval .and. mask_array(qdata%x,qdata%y) /= maskval) then
                  if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
                     distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) 
                     search_extrap = array(qdata%x, qdata%y, izz)
                  end if
               end if
                
            else
               if (array(qdata%x,qdata%y,izz) /= msgval) then
                  if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
                     distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) 
                     search_extrap = array(qdata%x, qdata%y, izz)
                  end if
               end if
            end if
         end do
      else
         search_extrap = msgval
      end if
  
      call q_destroy(q)
      call bitarray_destroy(b)
 
   end function search_extrap
 
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: four_pt_average
   !
   ! Purpose: Average of four surrounding grid point values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function four_pt_average(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to 
                                                !   interpolate within
      real, intent(in) :: xx, yy                ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: four_pt_average
  
      ! Local variables
      integer :: ifx, ify, icx, icy
      real :: fxfy, fxcy, cxfy, cxcy

      fxfy = 1.0
      fxcy = 1.0
      cxfy = 1.0
      cxcy = 1.0
  
      ifx = floor(xx)
      icx = ceiling(xx)
      ify = floor(yy)
      icy = ceiling(yy)

      ! First, make sure that the point is contained in the source array
      if (ifx < start_x .or. icx > end_x .or. &
          ify < start_y .or. icy > end_y) then

         ! But if the point is at most half a grid point out, we can
         !   still proceed with modified ifx, icx, ify, and icy.
         if (xx > real(start_x)-0.5 .and. ifx < start_x) then
            ifx = start_x
            icx = start_x
         else if (xx < real(end_x)+0.5 .and. icx > end_x) then
            ifx = end_x
            icx = end_x
         end if

         if (yy > real(start_y)-0.5 .and. ify < start_y) then
            ify = start_y
            icy = start_y
         else if (yy < real(end_y)+0.5 .and. icy > end_y) then
            ify = end_y
            icy = end_y
         end if

         if (ifx < start_x .or. icx > end_x .or. &
             ify < start_y .or. icy > end_y) then
            four_pt_average = msgval
            return
         end if
      end if

      if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
         ! we determine which maskval is useable by... if the symbol > is found, then only 
         !    values less than 'maskval' can be used and if the symbol < is found,  
         !    then only the values greater than 'maskval' can be used.
         if (mask_relational == '>') then
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0 
         else if (mask_relational == '<') then
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
         else  !equal
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0 
         end if 
      else if (present(mask_array) .and. present(maskval)) then
         if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0 
         if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0 
         if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0 
         if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0 
      else
         if (array(ifx, ify, izz) == msgval) fxfy = 0.0 
         if (array(ifx, icy, izz) == msgval) fxcy = 0.0 
         if (array(icx, ify, izz) == msgval) cxfy = 0.0 
         if (array(icx, icy, izz) == msgval) cxcy = 0.0 
      end if
  
      ! If all four points are missing, try the next interpolation method in the sequence
      if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
         four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                             start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
      else
         four_pt_average = (fxfy * array(ifx, ify, izz) + &
                            fxcy * array(ifx, icy, izz) + &
                            cxfy * array(icx, ify, izz) + &
                            cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
      end if
 
   end function four_pt_average


   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: wt_four_pt_average
   !
   ! Purpose: Weighted average of four surrounding grid point values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to 
                                                !   interpolate within
      real, intent(in) :: xx, yy                ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: wt_four_pt_average
  
      ! Local variables
      integer :: ifx, ify, icx, icy
      real :: fxfy, fxcy, cxfy, cxcy

      ifx = floor(xx)
      icx = ceiling(xx)
      ify = floor(yy)
      icy = ceiling(yy)

      fxfy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(ify))**2))
      fxcy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(icy))**2))
      cxfy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(ify))**2))
      cxcy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(icy))**2))

      ! First, make sure that the point is contained in the source array
      if (ifx < start_x .or. icx > end_x .or. &
          ify < start_y .or. icy > end_y) then

         ! But if the point is at most half a grid point out, we can
         !   still proceed with modified ifx, icx, ify, and icy.
         if (xx > real(start_x)-0.5 .and. ifx < start_x) then
            ifx = start_x
            icx = start_x
         else if (xx < real(end_x)+0.5 .and. icx > end_x) then
            ifx = end_x
            icx = end_x
         end if

         if (yy > real(start_y)-0.5 .and. ifx < start_y) then
            ify = start_y
            icy = start_y
         else if (yy < real(end_y)+0.5 .and. icy > end_y) then
            ify = end_y
            icy = end_y
         end if

         if (ifx < start_x .or. icx > end_x .or. &
             ify < start_y .or. icy > end_y) then
            wt_four_pt_average = msgval
            return
         end if
      end if

      if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then 
         ! we determine which maskval is useable by... if the symbol > is found, then only 
         !    values less than 'maskval' can be used and if the symbol < is found,  
         !    then only the values greater than 'maskval' can be used.
         if (mask_relational == '>') then
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0 
         else if (mask_relational == '<') then
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
         else  !equal
            if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0 
            if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0 
            if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0 
            if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0 
         end if              
      else if (present(mask_array) .and. present(maskval)) then
         if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0 
         if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0 
         if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0 
         if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0 
      else
         if (array(ifx, ify, izz) == msgval) fxfy = 0.0 
         if (array(ifx, icy, izz) == msgval) fxcy = 0.0 
         if (array(icx, ify, izz) == msgval) cxfy = 0.0 
         if (array(icx, icy, izz) == msgval) cxcy = 0.0 
      end if

      ! If all four points are missing, try the next interpolation method in the sequence
      if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
         wt_four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
      else
         wt_four_pt_average = (fxfy * array(ifx, ify, izz) + &
                               fxcy * array(ifx, icy, izz) + &
                               cxfy * array(icx, ify, izz) + &
                               cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
      end if
 
   end function wt_four_pt_average
 
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: sixteen_pt_average
   !
   ! Purpose: Average of sixteen surrounding grid point values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to
                                                !   interpolate within
      real, intent(in) :: xx , yy               ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: sixteen_pt_average
  
      ! Local variables
      integer :: i, j, ifx, ify
      real :: sum, sum_weight
      real, dimension(4,4) :: weights

      ifx = floor(xx)
      ify = floor(yy)
  
      ! First see whether the point is far enough within the array to 
      !   allow for a sixteen point average.
      if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
          ify < start_y+1 .or. ify > end_y-2) then
         sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
         return
      end if
  
      sum_weight = 0.0
      do i=1,4
         do j=1,4
            if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
               if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
                   weights(i,j) = 0.0
               else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
                   weights(i,j) = 0.0
               else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
                   weights(i,j) = 0.0
               else
                   weights(i,j) = 1.0
               end if
               if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
            else if (present(mask_array) .and. present(maskval)) then
               if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
                  weights(i,j) = 0.0
               else
                  weights(i,j) = 1.0
               end if
            else
               if (array(ifx+3-i, ify+3-j, izz) == msgval) then
                  weights(i,j) = 0.0
               else
                  weights(i,j) = 1.0
               end if
            end if
    
            sum_weight = sum_weight + weights(i,j)
   
         end do
      end do
  
      ! If all points are missing, try the next interpolation method in the sequence
      if (sum_weight == 0.0) then
         sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
      else
         sum = 0.0
         do i=1,4
            do j=1,4
               sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
            end do
         end do
         sixteen_pt_average = sum / sum_weight
      end if
 
   end function sixteen_pt_average
 
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: wt_sixteen_pt_average
   !
   ! Purpose: Weighted average of sixteen surrounding grid point values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
              start_y, end_y, start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to
                                                !   interpolate within
      real, intent(in) :: xx , yy               ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: wt_sixteen_pt_average
  
      ! Local variables
      integer :: i, j, ifx, ify
      real :: sum, sum_weight
      real, dimension(4,4) :: weights

      ifx = floor(xx)
      ify = floor(yy)
  
      ! First see whether the point is far enough within the array to 
      !   allow for a sixteen point average.
      if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
          ify < start_y+1 .or. ify > end_y-2) then
         wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                   start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
         return
      end if
  
      sum_weight = 0.0
      do i=1,4
         do j=1,4  
            if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
               if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
                   weights(i,j) = 0.0
               else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
                   weights(i,j) = 0.0
               else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
                   weights(i,j) = 0.0
               else
                   weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
               end if
               if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
            else if (present(mask_array) .and. present(maskval)) then
               if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
                  weights(i,j) = 0.0
               else
                  weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
               end if
            else
               if (array(ifx+3-i, ify+3-j, izz) == msgval) then
                  weights(i,j) = 0.0
               else
                  weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
               end if
            end if
    
            sum_weight = sum_weight + weights(i,j)
   
         end do
      end do
  
      ! If all points are missing, try the next interpolation method in the sequence
      if (sum_weight == 0.0) then
         wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                   start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
      else
         sum = 0.0
         do i=1,4
            do j=1,4
               sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
            end do
         end do
         wt_sixteen_pt_average = sum / sum_weight
      end if
 
   end function wt_sixteen_pt_average
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: four_pt
   !
   ! Purpose: Bilinear interpolation among four grid values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function four_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                       start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
 
      implicit none
  
      ! Arguments
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      integer, intent(in) :: izz                ! The z-index of the 2d-array to 
                                                !   interpolate within
      real, intent(in) :: xx , yy               ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: four_pt
  
      ! Local variables
      integer :: min_x, min_y, max_x, max_y

      min_x = floor(xx)
      min_y = floor(yy)
      max_x = ceiling(xx)
      max_y = ceiling(yy)
  
      if (min_x < start_x .or. max_x > end_x) then
         four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                   msgval, interp_list, idx, mask_relational, maskval, mask_array)
         return
      end if
  
      if (min_y < start_y .or. max_y > end_y) then
         four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                   msgval, interp_list, idx, mask_relational, maskval, mask_array)
         return
      end if
  
      ! If we have a missing value, try the next interpolation method in the sequence
      if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
         ! we determine which maskval is useable by... if the symbol > is found, then only 
         !    values less than 'maskval' can be used and if the symbol < is found,  
         !    then only the values greater than 'maskval' can be used.
         if (mask_relational == '>' ) then                 
            if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) > maskval .or. &
                array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) > maskval .or. &
                array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) > maskval .or. &
                array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) > maskval) then 
               four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                     msgval, interp_list, idx, mask_relational, maskval, mask_array)
               return
            end if
         else if (mask_relational == '<' ) then
            if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) < maskval .or. &
                array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) < maskval .or. &
                array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) < maskval .or. &
                array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) < maskval) then 
               four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                         msgval, interp_list, idx, mask_relational, maskval, mask_array)
               return
             end if           
         else if (mask_relational == ' ' ) then
            if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
                array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
                array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
                array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then 
               four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                         msgval, interp_list, idx, mask_relational, maskval, mask_array)
               return
            end if
         end if
      else if (present(mask_array) .and. present(maskval)) then
         if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
             array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
             array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
             array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then 
            four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                      msgval, interp_list, idx, mask_relational, maskval, mask_array)
         end if
      else
         if (array(min_x,min_y,izz) == msgval .or. &
             array(max_x,min_y,izz) == msgval .or. &
             array(min_x,max_y,izz) == msgval .or. &
             array(max_x,max_y,izz) == msgval ) then 
            four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                start_z, end_z, msgval, interp_list, idx)
            return
         end if
      end if
  
      if (min_x == max_x) then
         if (min_y == max_y) then
            four_pt = array(min_x,min_y,izz)
         else
            four_pt = array(min_x,min_y,izz)*(real(max_y)-yy) + &
                      array(min_x,max_y,izz)*(yy-real(min_y)) 
         end if
      else if (min_y == max_y) then
         if (min_x == max_x) then
            four_pt = array(min_x,min_y,izz)
         else
            four_pt = array(min_x,min_y,izz)*(real(max_x)-xx) + &
                      array(max_x,min_y,izz)*(xx-real(min_x)) 
         end if
      else
         four_pt = (yy - min_y) * (array(min_x,max_y,izz)*(real(max_x)-xx) + &
                                   array(max_x,max_y,izz)*(xx-real(min_x))) + &
                   (max_y - yy) * (array(min_x,min_y,izz)*(real(max_x)-xx) + &
                                   array(max_x,min_y,izz)*(xx-real(min_x)));
      end if
 
   end function four_pt
 
 
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! Name: sixteen_pt
   !
   ! Purpose: Overlapping parabolic interpolation among sixteen grid values
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   recursive function sixteen_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
                                 msgval, interp_list, idx, mask_relational, maskval, mask_array)
    
      implicit none
     
      ! Arguments
      integer, intent(in) :: izz    ! z-index of 2d-array to interpolate within
      integer, intent(in) :: start_x, start_y, start_z
      integer, intent(in) :: end_x, end_y, end_z
      real, intent(in) :: xx , yy              ! The location to interpolate to
      real, intent(in) :: msgval
      real, intent(in), optional :: maskval
      real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
          intent(in) :: array 
      integer, dimension(:), intent(in) :: interp_list
      integer, intent(in) :: idx
      real, dimension(start_x:end_x, start_y:end_y), &
          intent(in), optional :: mask_array 
      character (len=1), intent(in), optional :: mask_relational

      ! Return value
      real :: sixteen_pt
  
      ! Local variables
      integer :: n , i , j , k , kk , l , ll
      real :: x , y , a , b , c , d , e , f , g , h
      real, dimension(4,4) :: stl
      logical :: is_masked

      is_masked = .false.

      if (int(xx) < start_x .or. int(xx) > end_x .or. &
          int(yy) < start_y .or. int(yy) > end_y) then
         sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                      start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
         return
      end if
  
      sixteen_pt = 0.0
      n = 0
      i = int(xx + 0.00001)
      j = int(yy + 0.00001)
      x = xx - i
      y = yy - j

      if ( ( abs(x) > 0.0001 ) .or. ( abs(y) > 0.0001 ) ) then
  
         loop_1 : do k = 1,4
            kk = i + k - 2
            if ( kk < start_x) then
               kk = start_x
            else if ( kk > end_x) then
               kk = end_x
            end if
            loop_2 : do l = 1,4
               stl(k,l) = 0.
               ll = j + l - 2
               if ( ll < start_y ) then
                  ll = start_y
               else if ( ll > end_y) then
                  ll = end_y
               end if
               stl(k,l) = array(kk,ll,izz)
               n = n + 1
               if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
                  if (mask_relational == '>' .and. mask_array(kk,ll) > maskval) then
                     is_masked = .true.
                  else if (mask_relational == '<' .and. mask_array(kk,ll) < maskval) then
                     is_masked = .true.
                  else if (mask_relational == ' ' .and. mask_array(kk,ll) == maskval) then
                     is_masked = .true.
                  end if
               else if (present(mask_array) .and. present(maskval)) then
                  if (mask_array(kk,ll) == maskval) is_masked = .true.
               end if
               if ( stl(k,l) == 0. .and. msgval /= 0.) then
                  stl(k,l) = 1.E-20
               end if
            end do loop_2
         end do loop_1
   
         ! If we have a missing value, try the next interpolation method in the sequence
         if (present(mask_array) .and. present(maskval)) then
            do k=1,4
               do l=1,4
                  if (stl(k,l) == msgval .or. is_masked) then
                     sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
                                                  start_z, end_z, msgval, interp_list, idx, mask_relational, maskval, mask_array)
                     return
                  end if
               end do
            end do
         else
            do k=1,4
               do l=1,4
                  if (stl(k,l) == msgval) then
                     sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, …