!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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, end_y, &
start_z, end_z, msgval, interp_list, idx)
return
end if
end do
end do
end if
a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
sixteen_pt = oned(y,a,b,c,d)
if (n /= 16) then
e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
sixteen_pt = (sixteen_pt+oned(x,e,f,g,h)) * 0.5
end if
if (sixteen_pt == 1.E-20) sixteen_pt = 0.
else
if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
sixteen_pt = array(i,j,izz)
else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
sixteen_pt = array(i,j,izz)
else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
sixteen_pt = array(i,j,izz)
else
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)
end if
else if (present(mask_array) .and. present(maskval)) then
if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
sixteen_pt = array(i,j,izz)
else
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)
end if
else
if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. array(i,j,izz) /= msgval) then
sixteen_pt = array(i,j,izz)
else
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)
end if
end if
end if
end function sixteen_pt
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: oned
!
! Purpose: 1-dimensional overlapping parabolic interpolation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function oned(x,a,b,c,d)
implicit none
! Arguments
real, intent(in) :: x,a,b,c,d
! Return value
real :: oned
oned = 0.
if ( x == 0. ) then
oned = b
else if ( x == 1. ) then
oned = c
end if
if (b*c /= 0.) then
if ( a*d == 0. ) then
if ( ( a == 0 ) .and. ( d == 0 ) ) then
oned = b*(1.0-x)+c*x
else if ( a /= 0. ) then
oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
else if ( d /= 0. ) then
oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
end if
else
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)))
end if
end if
end function oned
end module interp_module