| blob | dd85a67a84334977e89ed9c63eeb5035283104d9 |
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! MODULE INTERP_MODULE
3 !
4 ! This module provides routines for interpolation.
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module interp_module
8 use bitarray_module
9 use misc_definitions_module
10 use module_debug
11 use queue_module
13 contains
15 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 ! Name: interp_array_from_string
17 !
18 ! Purpose:
19 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20 function interp_array_from_string(interp_string)
22 implicit none
24 ! Arguments
25 character (len=*), intent(in) :: interp_string
27 ! Local variables
28 integer :: j, p1, p2, iend, num_methods
30 ! Return value
31 integer, pointer, dimension(:) :: interp_array_from_string
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
40 allocate(interp_array_from_string(num_methods+1))
41 interp_array_from_string = 0
43 iend = len_trim(interp_string)
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) then
78 interp_array_from_string(j) = SEARCH
79 j = j + 1
80 else
81 if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
82 call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
83 end if
84 p1 = p2 + 1
85 p2 = index(interp_string(p1:iend),'+') + p1 - 1
86 end do
88 p2 = iend+1
89 if (p1 < iend) then
90 if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
91 len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
92 interp_array_from_string(j) = N_NEIGHBOR
93 j = j + 1
94 else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
95 len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
96 interp_array_from_string(j) = AVERAGE4
97 j = j + 1
98 else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
99 len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
100 interp_array_from_string(j) = AVERAGE16
101 j = j + 1
102 else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
103 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
104 interp_array_from_string(j) = W_AVERAGE4
105 j = j + 1
106 else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
107 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
108 interp_array_from_string(j) = W_AVERAGE16
109 j = j + 1
110 else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
111 len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
112 interp_array_from_string(j) = FOUR_POINT
113 j = j + 1
114 else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
115 len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
116 interp_array_from_string(j) = SIXTEEN_POINT
117 j = j + 1
118 else if (index(interp_string(p1:),'search') /= 0) then
119 interp_array_from_string(j) = SEARCH
120 j = j + 1
121 else
122 if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
123 call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
124 end if
125 end if
127 return
129 end function interp_array_from_string
132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
133 ! Name: interp_options_from_string
134 !
135 ! Purpose:
136 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
137 function interp_options_from_string(interp_string)
139 implicit none
141 ! Arguments
142 character (len=*), intent(in) :: interp_string
144 ! Local variables
145 integer :: j, p1, p2, iend, num_methods, istatus
147 ! Return value
148 integer, pointer, dimension(:) :: interp_options_from_string
150 ! Get an idea of how many interpolation methods are in the string
151 ! so we can allocate an appropriately sized array
152 num_methods = 1
153 do j=1,len_trim(interp_string)
154 if (interp_string(j:j) == '+') num_methods = num_methods + 1
155 end do
157 allocate(interp_options_from_string(num_methods+1))
158 interp_options_from_string(:) = 0
160 iend = len_trim(interp_string)
162 p1 = 1
163 p2 = index(interp_string(1:iend),'+')
164 j = 1
165 do while(p2 >= p1)
166 if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
167 len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
168 interp_options_from_string(j) = 0
169 j = j + 1
170 else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
171 len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
172 interp_options_from_string(j) = 0
173 j = j + 1
174 else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
175 len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
176 interp_options_from_string(j) = 0
177 j = j + 1
178 else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
179 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
180 interp_options_from_string(j) = 0
181 j = j + 1
182 else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
183 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
184 interp_options_from_string(j) = 0
185 j = j + 1
186 else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
187 len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
188 interp_options_from_string(j) = 0
189 j = j + 1
190 else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
191 len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
192 interp_options_from_string(j) = 0
193 j = j + 1
194 else if (index(interp_string(p1:p2-1),'search') /= 0) then
195 call get_search_depth(interp_string(p1:p2-1), interp_options_from_string(j), istatus)
196 j = j + 1
197 else
198 if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
199 call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
200 end if
201 p1 = p2 + 1
202 p2 = index(interp_string(p1:iend),'+') + p1 - 1
203 end do
205 p2 = iend+1
206 if (p1 < iend) then
207 if (index(interp_string(p1:p2-1),'nearest_neighbor') /= 0 .and. &
208 len_trim(interp_string(p1:p2-1)) == len_trim('nearest_neighbor')) then
209 interp_options_from_string(j) = 0
210 j = j + 1
211 else if (index(interp_string(p1:p2-1),'average_4pt') /= 0 .and. &
212 len_trim(interp_string(p1:p2-1)) == len_trim('average_4pt')) then
213 interp_options_from_string(j) = 0
214 j = j + 1
215 else if (index(interp_string(p1:p2-1),'average_16pt') /= 0 .and. &
216 len_trim(interp_string(p1:p2-1)) == len_trim('average_16pt')) then
217 interp_options_from_string(j) = 0
218 j = j + 1
219 else if (index(interp_string(p1:p2-1),'wt_average_4pt') /= 0 .and. &
220 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_4pt')) then
221 interp_options_from_string(j) = 0
222 j = j + 1
223 else if (index(interp_string(p1:p2-1),'wt_average_16pt') /= 0 .and. &
224 len_trim(interp_string(p1:p2-1)) == len_trim('wt_average_16pt')) then
225 interp_options_from_string(j) = 0
226 j = j + 1
227 else if (index(interp_string(p1:p2-1),'four_pt') /= 0 .and. &
228 len_trim(interp_string(p1:p2-1)) == len_trim('four_pt')) then
229 interp_options_from_string(j) = 0
230 j = j + 1
231 else if (index(interp_string(p1:p2-1),'sixteen_pt') /= 0 .and. &
232 len_trim(interp_string(p1:p2-1)) == len_trim('sixteen_pt')) then
233 interp_options_from_string(j) = 0
234 j = j + 1
235 else if (index(interp_string(p1:),'search') /= 0) then
236 call get_search_depth(interp_string(p1:), interp_options_from_string(j), istatus)
237 j = j + 1
238 else
239 if (index(interp_string(p1:p2-1),'average_gcell') == 0) &
240 call mprintf(.true.,WARN,'Unrecognized interpolation method %s.',s1=interp_string(p1:p2-1))
241 end if
242 end if
244 return
246 end function interp_options_from_string
249 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
250 ! Name: get_search_depth
251 !
252 ! Pupose:
253 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254 subroutine get_search_depth(interp_opt, depth, istatus)
256 implicit none
258 ! Arguments
259 integer, intent(out) :: istatus
260 integer, intent(out) :: depth
261 character (len=*), intent(in) :: interp_opt
263 ! Local variables
264 integer :: i, p1, p2, p
266 istatus = 1
267 depth = 1200
269 i = index(interp_opt,'search')
270 if (i /= 0) then
272 ! Check for a max search depth
273 p = index(interp_opt(i:),'+')
274 if (p == 0) p = len_trim(interp_opt)
275 p1 = index(interp_opt(i:p),'(')
276 p2 = index(interp_opt(i:p),')')
277 if (p1 /= 0 .and. p2 /= 0) then
278 read(interp_opt(p1+1:p2-1),*,err=1000) depth
279 else if (p1 == 0 .and. p2 == 0) then
280 ! keep depth at 1200, no warning
281 else
282 call mprintf(.true., WARN, 'Problem with specified search depth '// &
283 'for search interp option. Setting max depth to 1200.')
284 depth = 1200
285 end if
286 end if
287 istatus = 0
289 return
291 1000 call mprintf(.true., ERROR, 'Search depth option to search interpolator '// &
292 'must be an integer value, enclosed in parentheses immediately '// &
293 'after keyword "search"')
295 end subroutine get_search_depth
298 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
299 ! Name: interp_sequence
300 !
301 ! Purpose: Delegates the actual task of interpolation to specific
302 ! interpolation routines defined in the module.
303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
304 recursive function interp_sequence(xx, yy, izz, array, start_x, end_x, &
305 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
307 implicit none
309 ! Arguments
310 integer, intent(in) :: start_x, start_y, start_z
311 integer, intent(in) :: end_x, end_y, end_z
312 integer, intent(in) :: izz ! The z-index of the 2d-array to
313 ! interpolate within
314 real, intent(in) :: xx , yy ! The location to interpolate to
315 real, intent(in) :: msgval
316 real, intent(in), optional :: maskval
317 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
318 integer, intent(in) :: idx
319 integer, dimension(:), intent(in) :: interp_list
320 integer, dimension(:), intent(in) :: interp_opts
321 real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
322 character (len=1), intent(in), optional :: mask_relational
324 ! Return value
325 real :: interp_sequence
327 ! No more interpolation methods to try
328 if (interp_list(idx) == 0) then
329 interp_sequence = msgval
330 return
331 end if
333 if (interp_list(idx) == FOUR_POINT) then
334 interp_sequence = four_pt(xx, yy, izz, array, start_x, end_x, &
335 start_y, end_y, start_z, end_z, &
336 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
337 else if (interp_list(idx) == AVERAGE4) then
338 interp_sequence = four_pt_average(xx, yy, izz, array, start_x, end_x, &
339 start_y, end_y, start_z, end_z, &
340 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
341 else if (interp_list(idx) == W_AVERAGE4) then
342 interp_sequence = wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
343 start_y, end_y, start_z, end_z, &
344 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
345 else if (interp_list(idx) == N_NEIGHBOR) then
346 interp_sequence = nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
347 start_y, end_y, start_z, end_z, &
348 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
349 else if (interp_list(idx) == SIXTEEN_POINT) then
350 interp_sequence = sixteen_pt(xx, yy, izz, array, start_x, end_x, &
351 start_y, end_y, start_z, end_z, &
352 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
353 else if (interp_list(idx) == SEARCH) then
354 interp_sequence = search_extrap(xx, yy, izz, array, start_x, end_x, &
355 start_y, end_y, start_z, end_z, &
356 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
357 else if (interp_list(idx) == AVERAGE16) then
358 interp_sequence = sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
359 start_y, end_y, start_z, end_z, &
360 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
361 else if (interp_list(idx) == W_AVERAGE16) then
362 interp_sequence = wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
363 start_y, end_y, start_z, end_z, &
364 msgval, interp_list, interp_opts, idx+1, mask_relational, maskval, mask_array)
365 end if
367 end function interp_sequence
370 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
371 ! Name: nearest_neighbor
372 !
373 ! Purpose: Returns the point nearest to (xx,yy). If (xx,yy) is outside of the
374 ! array, the point on the edge of the array nearest to (xx,yy) is returned.
375 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
376 recursive function nearest_neighbor(xx, yy, izz, array, start_x, end_x, &
377 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
379 implicit none
381 ! Arguments
382 integer, intent(in) :: start_x, start_y, start_z
383 integer, intent(in) :: end_x, end_y, end_z
384 integer, intent(in) :: izz ! The z-index of the 2d-array to
385 ! interpolate within
386 real, intent(in) :: xx , yy ! The location to interpolate to
387 real, intent(in) :: msgval
388 real, intent(in), optional :: maskval
389 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
390 integer, dimension(:), intent(in) :: interp_list
391 integer, dimension(:), intent(in) :: interp_opts
392 integer, intent(in) :: idx
393 real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
394 character (len=1), intent(in), optional :: mask_relational
396 ! Return value
397 real :: nearest_neighbor
399 ! Local variables
400 integer :: ix, iy
402 ix = nint(xx)
403 iy = nint(yy)
405 ! The first thing to do is to ensure that the point (xx,yy) is within the array
406 if (ix < start_x .or. ix > end_x) then
407 nearest_neighbor = msgval
408 return
409 end if
411 if (iy < start_y .or. iy > end_y) then
412 nearest_neighbor = msgval
413 return
414 end if
416 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
417 if (mask_relational == '<' .and. mask_array(ix,iy) < maskval) then
418 nearest_neighbor = msgval
419 else if (mask_relational == '>' .and. mask_array(ix,iy) > maskval) then
420 nearest_neighbor = msgval
421 else if (mask_relational == ' ' .and. mask_array(ix,iy) == maskval) then
422 nearest_neighbor = msgval
423 else
424 nearest_neighbor = array(ix,iy,izz)
425 end if
426 else if (present(mask_array) .and. present(maskval)) then
427 if (maskval == mask_array(ix,iy)) then
428 nearest_neighbor = msgval
429 else
430 nearest_neighbor = array(ix,iy,izz)
431 end if
432 else
433 nearest_neighbor = array(ix,iy,izz)
434 end if
436 ! If we have a missing value, try the next interpolation method in the sequence
437 if (nearest_neighbor == msgval) then
438 nearest_neighbor = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
439 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
440 end if
442 end function nearest_neighbor
445 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
446 ! Name: search_extrap
447 !
448 ! Purpose: Returns the point nearest to (xx,yy) that has a non-missing value.
449 ! If no valid value can be found in the array, msgval is returned.
450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
451 recursive function search_extrap(xx, yy, izz, array, start_x, end_x, &
452 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
454 implicit none
456 ! Arguments
457 integer, intent(in) :: start_x, start_y, start_z
458 integer, intent(in) :: end_x, end_y, end_z
459 integer, intent(in) :: izz ! The z-index of the 2d-array to search
460 real, intent(in) :: xx , yy ! The location of the search origin
461 real, intent(in) :: msgval
462 real, intent(in), optional :: maskval
463 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
464 intent(in) :: array
465 integer, dimension(:), intent(in) :: interp_list
466 integer, dimension(:), intent(in) :: interp_opts
467 integer, intent(in) :: idx
468 real, dimension(start_x:end_x, start_y:end_y), &
469 intent(in), optional :: mask_array
470 character (len=1), intent(in), optional :: mask_relational
472 ! Return value
473 real :: search_extrap
475 ! Local variables
476 integer :: i, j
477 real :: distance
478 logical :: found_valid
479 type (bitarray) :: b
480 type (queue) :: q
481 type (q_data) :: qdata
483 ! We only search if the starting point is within the array
484 if (nint(xx) < start_x .or. nint(xx) > end_x .or. &
485 nint(yy) < start_y .or. nint(yy) > end_y) then
486 search_extrap = msgval
487 return
488 end if
490 call bitarray_create(b, (end_x-start_x+1), (end_y-start_y+1))
491 call q_init(q)
493 found_valid = .false.
494 qdata%x = nint(xx)
495 qdata%y = nint(yy)
496 qdata%depth = 0
497 call q_insert(q, qdata)
498 call bitarray_set(b, qdata%x-start_x+1, qdata%y-start_y+1)
500 do while (q_isdata(q) .and. (.not. found_valid))
501 qdata = q_remove(q)
502 i = qdata%x
503 j = qdata%y
505 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
506 if (mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
507 found_valid = .true.
508 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
509 found_valid = .true.
510 else if (mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
511 found_valid = .true.
512 end if
513 else if (present(mask_array) .and. present(maskval)) then
514 if (array(i,j,izz) /= msgval .and. mask_array(i,j) /= maskval) found_valid = .true.
515 else
516 if (array(i,j,izz) /= msgval) found_valid = .true.
517 end if
519 if (i-1 >= start_x) then
520 if (.not. bitarray_test(b, (i-1)-start_x+1, j-start_y+1)) then
521 if (qdata%depth < interp_opts(idx-1)) then ! idx-1, since idx was incremented before call to this subroutine
522 qdata%x = i-1
523 qdata%y = j
524 qdata%depth = qdata%depth+1
525 call q_insert(q, qdata)
526 call bitarray_set(b, (i-1)-start_x+1, j-start_y+1)
527 end if
528 end if
529 end if
530 if (i+1 <= end_x) then
531 if (.not. bitarray_test(b, (i+1)-start_x+1, j-start_y+1)) then
532 if (qdata%depth < interp_opts(idx-1)) then ! idx-1, since idx was incremented before call to this subroutine
533 qdata%x = i+1
534 qdata%y = j
535 qdata%depth = qdata%depth+1
536 call q_insert(q, qdata)
537 call bitarray_set(b, (i+1)-start_x+1, j-start_y+1)
538 end if
539 end if
540 end if
541 if (j-1 >= start_y) then
542 if (.not. bitarray_test(b, i-start_x+1, (j-1)-start_y+1)) then
543 if (qdata%depth < interp_opts(idx-1)) then ! idx-1, since idx was incremented before call to this subroutine
544 qdata%x = i
545 qdata%y = j-1
546 qdata%depth = qdata%depth+1
547 call q_insert(q, qdata)
548 call bitarray_set(b, i-start_x+1, (j-1)-start_y+1)
549 end if
550 end if
551 end if
552 if (j+1 <= end_y) then
553 if (.not. bitarray_test(b, i-start_x+1, (j+1)-start_y+1)) then
554 if (qdata%depth < interp_opts(idx-1)) then ! idx-1, since idx was incremented before call to this subroutine
555 qdata%x = i
556 qdata%y = j+1
557 qdata%depth = qdata%depth+1
558 call q_insert(q, qdata)
559 call bitarray_set(b, i-start_x+1, (j+1)-start_y+1)
560 end if
561 end if
562 end if
563 end do
565 if (found_valid) then
566 distance = (real(i)-xx)*(real(i)-xx)+(real(j)-yy)*(real(j)-yy)
567 search_extrap = array(i,j,izz)
568 do while (q_isdata(q))
569 qdata = q_remove(q)
570 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
571 if (mask_relational == '<' .and. mask_array(qdata%x,qdata%y) >= maskval &
572 .and. array(qdata%x,qdata%y,izz) /= msgval) then
573 if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
574 distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
575 search_extrap = array(qdata%x, qdata%y, izz)
576 end if
577 else if (mask_relational == '>' .and. mask_array(qdata%x,qdata%y) <= maskval &
578 .and. array(qdata%x,qdata%y,izz) /= msgval) then
579 if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
580 distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
581 search_extrap = array(qdata%x, qdata%y, izz)
582 end if
583 else if (mask_relational == ' ' .and. mask_array(qdata%x,qdata%y) /= maskval &
584 .and. array(qdata%x,qdata%y,izz) /= msgval) then
585 if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
586 distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
587 search_extrap = array(qdata%x, qdata%y, izz)
588 end if
589 end if
591 else if (present(mask_array) .and. present(maskval)) then
592 if (array(qdata%x,qdata%y,izz) /= msgval .and. mask_array(qdata%x,qdata%y) /= maskval) then
593 if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
594 distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
595 search_extrap = array(qdata%x, qdata%y, izz)
596 end if
597 end if
599 else
600 if (array(qdata%x,qdata%y,izz) /= msgval) then
601 if ((real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy) < distance) then
602 distance = (real(qdata%x)-xx)*(real(qdata%x)-xx)+(real(qdata%y)-yy)*(real(qdata%y)-yy)
603 search_extrap = array(qdata%x, qdata%y, izz)
604 end if
605 end if
606 end if
607 end do
608 else
609 search_extrap = msgval
610 end if
612 call q_destroy(q)
613 call bitarray_destroy(b)
615 end function search_extrap
618 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
619 ! Name: four_pt_average
620 !
621 ! Purpose: Average of four surrounding grid point values
622 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
623 recursive function four_pt_average(xx, yy, izz, array, start_x, end_x, &
624 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
626 implicit none
628 ! Arguments
629 integer, intent(in) :: start_x, start_y, start_z
630 integer, intent(in) :: end_x, end_y, end_z
631 integer, intent(in) :: izz ! The z-index of the 2d-array to
632 ! interpolate within
633 real, intent(in) :: xx, yy ! The location to interpolate to
634 real, intent(in) :: msgval
635 real, intent(in), optional :: maskval
636 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
637 intent(in) :: array
638 integer, dimension(:), intent(in) :: interp_list
639 integer, dimension(:), intent(in) :: interp_opts
640 integer, intent(in) :: idx
641 real, dimension(start_x:end_x, start_y:end_y), &
642 intent(in), optional :: mask_array
643 character (len=1), intent(in), optional :: mask_relational
645 ! Return value
646 real :: four_pt_average
648 ! Local variables
649 integer :: ifx, ify, icx, icy
650 real :: fxfy, fxcy, cxfy, cxcy
652 fxfy = 1.0
653 fxcy = 1.0
654 cxfy = 1.0
655 cxcy = 1.0
657 ifx = floor(xx)
658 icx = ceiling(xx)
659 ify = floor(yy)
660 icy = ceiling(yy)
662 ! First, make sure that the point is contained in the source array
663 if (ifx < start_x .or. icx > end_x .or. &
664 ify < start_y .or. icy > end_y) then
666 ! But if the point is at most half a grid point out, we can
667 ! still proceed with modified ifx, icx, ify, and icy.
668 if (xx > real(start_x)-0.5 .and. ifx < start_x) then
669 ifx = start_x
670 icx = start_x
671 else if (xx < real(end_x)+0.5 .and. icx > end_x) then
672 ifx = end_x
673 icx = end_x
674 end if
676 if (yy > real(start_y)-0.5 .and. ify < start_y) then
677 ify = start_y
678 icy = start_y
679 else if (yy < real(end_y)+0.5 .and. icy > end_y) then
680 ify = end_y
681 icy = end_y
682 end if
684 if (ifx < start_x .or. icx > end_x .or. &
685 ify < start_y .or. icy > end_y) then
686 four_pt_average = msgval
687 return
688 end if
689 end if
691 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
692 ! we determine which maskval is useable by... if the symbol > is found, then only
693 ! values less than 'maskval' can be used and if the symbol < is found,
694 ! then only the values greater than 'maskval' can be used.
695 if (mask_relational == '>') then
696 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0
697 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0
698 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0
699 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0
700 else if (mask_relational == '<') then
701 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0
702 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0
703 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0
704 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
705 else !equal
706 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
707 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
708 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
709 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
710 end if
711 else if (present(mask_array) .and. present(maskval)) then
712 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
713 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
714 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
715 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
716 else
717 if (array(ifx, ify, izz) == msgval) fxfy = 0.0
718 if (array(ifx, icy, izz) == msgval) fxcy = 0.0
719 if (array(icx, ify, izz) == msgval) cxfy = 0.0
720 if (array(icx, icy, izz) == msgval) cxcy = 0.0
721 end if
723 ! If all four points are missing, try the next interpolation method in the sequence
724 if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
725 four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
726 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
727 else
728 four_pt_average = (fxfy * array(ifx, ify, izz) + &
729 fxcy * array(ifx, icy, izz) + &
730 cxfy * array(icx, ify, izz) + &
731 cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
732 end if
734 end function four_pt_average
737 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
738 ! Name: wt_four_pt_average
739 !
740 ! Purpose: Weighted average of four surrounding grid point values
741 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
742 recursive function wt_four_pt_average(xx, yy, izz, array, start_x, end_x, &
743 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
745 implicit none
747 ! Arguments
748 integer, intent(in) :: start_x, start_y, start_z
749 integer, intent(in) :: end_x, end_y, end_z
750 integer, intent(in) :: izz ! The z-index of the 2d-array to
751 ! interpolate within
752 real, intent(in) :: xx, yy ! The location to interpolate to
753 real, intent(in) :: msgval
754 real, intent(in), optional :: maskval
755 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
756 intent(in) :: array
757 integer, dimension(:), intent(in) :: interp_list
758 integer, dimension(:), intent(in) :: interp_opts
759 integer, intent(in) :: idx
760 real, dimension(start_x:end_x, start_y:end_y), &
761 intent(in), optional :: mask_array
762 character (len=1), intent(in), optional :: mask_relational
764 ! Return value
765 real :: wt_four_pt_average
767 ! Local variables
768 integer :: ifx, ify, icx, icy
769 real :: fxfy, fxcy, cxfy, cxcy
771 ifx = floor(xx)
772 icx = ceiling(xx)
773 ify = floor(yy)
774 icy = ceiling(yy)
776 fxfy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(ify))**2))
777 fxcy = max(0., 1.0 - sqrt((xx-real(ifx))**2+(yy-real(icy))**2))
778 cxfy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(ify))**2))
779 cxcy = max(0., 1.0 - sqrt((xx-real(icx))**2+(yy-real(icy))**2))
781 ! First, make sure that the point is contained in the source array
782 if (ifx < start_x .or. icx > end_x .or. &
783 ify < start_y .or. icy > end_y) then
785 ! But if the point is at most half a grid point out, we can
786 ! still proceed with modified ifx, icx, ify, and icy.
787 if (xx > real(start_x)-0.5 .and. ifx < start_x) then
788 ifx = start_x
789 icx = start_x
790 else if (xx < real(end_x)+0.5 .and. icx > end_x) then
791 ifx = end_x
792 icx = end_x
793 end if
795 if (yy > real(start_y)-0.5 .and. ifx < start_y) then
796 ify = start_y
797 icy = start_y
798 else if (yy < real(end_y)+0.5 .and. icy > end_y) then
799 ify = end_y
800 icy = end_y
801 end if
803 if (ifx < start_x .or. icx > end_x .or. &
804 ify < start_y .or. icy > end_y) then
805 wt_four_pt_average = msgval
806 return
807 end if
808 end if
810 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
811 ! we determine which maskval is useable by... if the symbol > is found, then only
812 ! values less than 'maskval' can be used and if the symbol < is found,
813 ! then only the values greater than 'maskval' can be used.
814 if (mask_relational == '>') then
815 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) > maskval) fxfy = 0.0
816 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) > maskval) fxcy = 0.0
817 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) > maskval) cxfy = 0.0
818 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) > maskval) cxcy = 0.0
819 else if (mask_relational == '<') then
820 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) < maskval) fxfy = 0.0
821 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) < maskval) fxcy = 0.0
822 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) < maskval) cxfy = 0.0
823 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) < maskval) cxcy = 0.0
824 else !equal
825 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
826 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
827 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
828 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
829 end if
830 else if (present(mask_array) .and. present(maskval)) then
831 if (array(ifx, ify, izz) == msgval .or. mask_array(ifx,ify) == maskval) fxfy = 0.0
832 if (array(ifx, icy, izz) == msgval .or. mask_array(ifx,icy) == maskval) fxcy = 0.0
833 if (array(icx, ify, izz) == msgval .or. mask_array(icx,ify) == maskval) cxfy = 0.0
834 if (array(icx, icy, izz) == msgval .or. mask_array(icx,icy) == maskval) cxcy = 0.0
835 else
836 if (array(ifx, ify, izz) == msgval) fxfy = 0.0
837 if (array(ifx, icy, izz) == msgval) fxcy = 0.0
838 if (array(icx, ify, izz) == msgval) cxfy = 0.0
839 if (array(icx, icy, izz) == msgval) cxcy = 0.0
840 end if
842 ! If all four points are missing, try the next interpolation method in the sequence
843 if (fxfy == 0.0 .and. fxcy == 0.0 .and. cxfy == 0.0 .and. cxcy == 0.0) then
844 wt_four_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
845 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
846 else
847 wt_four_pt_average = (fxfy * array(ifx, ify, izz) + &
848 fxcy * array(ifx, icy, izz) + &
849 cxfy * array(icx, ify, izz) + &
850 cxcy * array(icx, icy, izz) ) / (fxfy + fxcy + cxfy + cxcy)
851 end if
853 end function wt_four_pt_average
856 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
857 ! Name: sixteen_pt_average
858 !
859 ! Purpose: Average of sixteen surrounding grid point values
860 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
861 recursive function sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
862 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
864 implicit none
866 ! Arguments
867 integer, intent(in) :: start_x, start_y, start_z
868 integer, intent(in) :: end_x, end_y, end_z
869 integer, intent(in) :: izz ! The z-index of the 2d-array to
870 ! interpolate within
871 real, intent(in) :: xx , yy ! The location to interpolate to
872 real, intent(in) :: msgval
873 real, intent(in), optional :: maskval
874 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
875 intent(in) :: array
876 integer, dimension(:), intent(in) :: interp_list
877 integer, dimension(:), intent(in) :: interp_opts
878 integer, intent(in) :: idx
879 real, dimension(start_x:end_x, start_y:end_y), &
880 intent(in), optional :: mask_array
881 character (len=1), intent(in), optional :: mask_relational
883 ! Return value
884 real :: sixteen_pt_average
886 ! Local variables
887 integer :: i, j, ifx, ify
888 real :: sum, sum_weight
889 real, dimension(4,4) :: weights
891 ifx = floor(xx)
892 ify = floor(yy)
894 ! First see whether the point is far enough within the array to
895 ! allow for a sixteen point average.
896 if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
897 ify < start_y+1 .or. ify > end_y-2) then
898 sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
899 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
900 return
901 end if
903 sum_weight = 0.0
904 do i=1,4
905 do j=1,4
906 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
907 if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
908 weights(i,j) = 0.0
909 else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
910 weights(i,j) = 0.0
911 else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
912 weights(i,j) = 0.0
913 else
914 weights(i,j) = 1.0
915 end if
916 if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
917 else if (present(mask_array) .and. present(maskval)) then
918 if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
919 weights(i,j) = 0.0
920 else
921 weights(i,j) = 1.0
922 end if
923 else
924 if (array(ifx+3-i, ify+3-j, izz) == msgval) then
925 weights(i,j) = 0.0
926 else
927 weights(i,j) = 1.0
928 end if
929 end if
931 sum_weight = sum_weight + weights(i,j)
933 end do
934 end do
936 ! If all points are missing, try the next interpolation method in the sequence
937 if (sum_weight == 0.0) then
938 sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
939 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
940 else
941 sum = 0.0
942 do i=1,4
943 do j=1,4
944 sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
945 end do
946 end do
947 sixteen_pt_average = sum / sum_weight
948 end if
950 end function sixteen_pt_average
953 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
954 ! Name: wt_sixteen_pt_average
955 !
956 ! Purpose: Weighted average of sixteen surrounding grid point values
957 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
958 recursive function wt_sixteen_pt_average(xx, yy, izz, array, start_x, end_x, &
959 start_y, end_y, start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
961 implicit none
963 ! Arguments
964 integer, intent(in) :: start_x, start_y, start_z
965 integer, intent(in) :: end_x, end_y, end_z
966 integer, intent(in) :: izz ! The z-index of the 2d-array to
967 ! interpolate within
968 real, intent(in) :: xx , yy ! The location to interpolate to
969 real, intent(in) :: msgval
970 real, intent(in), optional :: maskval
971 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
972 intent(in) :: array
973 integer, dimension(:), intent(in) :: interp_list
974 integer, dimension(:), intent(in) :: interp_opts
975 integer, intent(in) :: idx
976 real, dimension(start_x:end_x, start_y:end_y), &
977 intent(in), optional :: mask_array
978 character (len=1), intent(in), optional :: mask_relational
980 ! Return value
981 real :: wt_sixteen_pt_average
983 ! Local variables
984 integer :: i, j, ifx, ify
985 real :: sum, sum_weight
986 real, dimension(4,4) :: weights
988 ifx = floor(xx)
989 ify = floor(yy)
991 ! First see whether the point is far enough within the array to
992 ! allow for a sixteen point average.
993 if (ifx < start_x+1 .or. ifx > end_x-2 .or. &
994 ify < start_y+1 .or. ify > end_y-2) then
995 wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
996 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
997 return
998 end if
1000 sum_weight = 0.0
1001 do i=1,4
1002 do j=1,4
1003 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1004 if (mask_relational == '>' .and. mask_array(ifx+3-i, ify+3-j) > maskval) then
1005 weights(i,j) = 0.0
1006 else if (mask_relational == '<' .and. mask_array(ifx+3-i, ify+3-j) < maskval) then
1007 weights(i,j) = 0.0
1008 else if (mask_relational == ' ' .and. mask_array(ifx+3-i, ify+3-j) == maskval) then
1009 weights(i,j) = 0.0
1010 else
1011 weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
1012 end if
1013 if (array(ifx+3-i, ify+3-j, izz) == msgval) weights(i,j) = 0.0
1014 else if (present(mask_array) .and. present(maskval)) then
1015 if (array(ifx+3-i, ify+3-j, izz) == msgval .or. mask_array(ifx+3-i, ify+3-j) == maskval) then
1016 weights(i,j) = 0.0
1017 else
1018 weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
1019 end if
1020 else
1021 if (array(ifx+3-i, ify+3-j, izz) == msgval) then
1022 weights(i,j) = 0.0
1023 else
1024 weights(i,j) = max(0., 2.0 - sqrt((xx-real(ifx+3-i))**2+(yy-real(ify+3-j))**2))
1025 end if
1026 end if
1028 sum_weight = sum_weight + weights(i,j)
1030 end do
1031 end do
1033 ! If all points are missing, try the next interpolation method in the sequence
1034 if (sum_weight == 0.0) then
1035 wt_sixteen_pt_average = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1036 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1037 else
1038 sum = 0.0
1039 do i=1,4
1040 do j=1,4
1041 sum = sum + weights(i,j) * array(ifx+3-i, ify+3-j, izz)
1042 end do
1043 end do
1044 wt_sixteen_pt_average = sum / sum_weight
1045 end if
1047 end function wt_sixteen_pt_average
1050 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1051 ! Name: four_pt
1052 !
1053 ! Purpose: Bilinear interpolation among four grid values
1054 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1055 recursive function four_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1056 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1058 implicit none
1060 ! Arguments
1061 integer, intent(in) :: start_x, start_y, start_z
1062 integer, intent(in) :: end_x, end_y, end_z
1063 integer, intent(in) :: izz ! The z-index of the 2d-array to
1064 ! interpolate within
1065 real, intent(in) :: xx , yy ! The location to interpolate to
1066 real, intent(in) :: msgval
1067 real, intent(in), optional :: maskval
1068 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), intent(in) :: array
1069 integer, dimension(:), intent(in) :: interp_list
1070 integer, dimension(:), intent(in) :: interp_opts
1071 integer, intent(in) :: idx
1072 real, dimension(start_x:end_x, start_y:end_y), intent(in), optional :: mask_array
1073 character (len=1), intent(in), optional :: mask_relational
1075 ! Return value
1076 real :: four_pt
1078 ! Local variables
1079 integer :: min_x, min_y, max_x, max_y
1081 min_x = floor(xx)
1082 min_y = floor(yy)
1083 max_x = ceiling(xx)
1084 max_y = ceiling(yy)
1086 if (min_x < start_x .or. max_x > end_x) then
1087 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1088 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1089 return
1090 end if
1092 if (min_y < start_y .or. max_y > end_y) then
1093 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1094 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1095 return
1096 end if
1098 ! If we have a missing value, try the next interpolation method in the sequence
1099 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1100 ! we determine which maskval is useable by... if the symbol > is found, then only
1101 ! values less than 'maskval' can be used and if the symbol < is found,
1102 ! then only the values greater than 'maskval' can be used.
1103 if (mask_relational == '>' ) then
1104 if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) > maskval .or. &
1105 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) > maskval .or. &
1106 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) > maskval .or. &
1107 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) > maskval) then
1108 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1109 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1110 return
1111 end if
1112 else if (mask_relational == '<' ) then
1113 if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) < maskval .or. &
1114 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) < maskval .or. &
1115 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) < maskval .or. &
1116 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) < maskval) then
1117 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1118 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1119 return
1120 end if
1121 else if (mask_relational == ' ' ) then
1122 if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
1123 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
1124 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
1125 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then
1126 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1127 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1128 return
1129 end if
1130 end if
1131 else if (present(mask_array) .and. present(maskval)) then
1132 if (array(min_x,min_y,izz) == msgval .or. mask_array(min_x,min_y) == maskval .or. &
1133 array(max_x,min_y,izz) == msgval .or. mask_array(max_x,min_y) == maskval .or. &
1134 array(min_x,max_y,izz) == msgval .or. mask_array(min_x,max_y) == maskval .or. &
1135 array(max_x,max_y,izz) == msgval .or. mask_array(max_x,max_y) == maskval) then
1136 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1137 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1138 end if
1139 else
1140 if (array(min_x,min_y,izz) == msgval .or. &
1141 array(max_x,min_y,izz) == msgval .or. &
1142 array(min_x,max_y,izz) == msgval .or. &
1143 array(max_x,max_y,izz) == msgval ) then
1144 four_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1145 start_z, end_z, msgval, interp_list, interp_opts, idx)
1146 return
1147 end if
1148 end if
1150 if (min_x == max_x) then
1151 if (min_y == max_y) then
1152 four_pt = array(min_x,min_y,izz)
1153 else
1154 four_pt = array(min_x,min_y,izz)*(real(max_y)-yy) + &
1155 array(min_x,max_y,izz)*(yy-real(min_y))
1156 end if
1157 else if (min_y == max_y) then
1158 if (min_x == max_x) then
1159 four_pt = array(min_x,min_y,izz)
1160 else
1161 four_pt = array(min_x,min_y,izz)*(real(max_x)-xx) + &
1162 array(max_x,min_y,izz)*(xx-real(min_x))
1163 end if
1164 else
1165 four_pt = (yy - min_y) * (array(min_x,max_y,izz)*(real(max_x)-xx) + &
1166 array(max_x,max_y,izz)*(xx-real(min_x))) + &
1167 (max_y - yy) * (array(min_x,min_y,izz)*(real(max_x)-xx) + &
1168 array(max_x,min_y,izz)*(xx-real(min_x)));
1169 end if
1171 end function four_pt
1174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1175 ! Name: sixteen_pt
1176 !
1177 ! Purpose: Overlapping parabolic interpolation among sixteen grid values
1178 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 recursive function sixteen_pt(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1180 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1182 implicit none
1184 ! Arguments
1185 integer, intent(in) :: izz ! z-index of 2d-array to interpolate within
1186 integer, intent(in) :: start_x, start_y, start_z
1187 integer, intent(in) :: end_x, end_y, end_z
1188 real, intent(in) :: xx , yy ! The location to interpolate to
1189 real, intent(in) :: msgval
1190 real, intent(in), optional :: maskval
1191 real, dimension(start_x:end_x, start_y:end_y, start_z:end_z), &
1192 intent(in) :: array
1193 integer, dimension(:), intent(in) :: interp_list
1194 integer, dimension(:), intent(in) :: interp_opts
1195 integer, intent(in) :: idx
1196 real, dimension(start_x:end_x, start_y:end_y), &
1197 intent(in), optional :: mask_array
1198 character (len=1), intent(in), optional :: mask_relational
1200 ! Return value
1201 real :: sixteen_pt
1203 ! Local variables
1204 integer :: n , i , j , k , kk , l , ll
1205 real :: x , y , a , b , c , d , e , f , g , h
1206 real, dimension(4,4) :: stl
1207 logical :: is_masked
1209 is_masked = .false.
1211 if (int(xx) < start_x .or. int(xx) > end_x .or. &
1212 int(yy) < start_y .or. int(yy) > end_y) then
1213 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1214 start_z, end_z, msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1215 return
1216 end if
1218 sixteen_pt = 0.0
1219 n = 0
1220 i = int(xx + 0.00001)
1221 j = int(yy + 0.00001)
1222 x = xx - i
1223 y = yy - j
1225 if ( ( abs(x) > 0.0001 ) .or. ( abs(y) > 0.0001 ) ) then
1227 loop_1 : do k = 1,4
1228 kk = i + k - 2
1229 if ( kk < start_x) then
1230 kk = start_x
1231 else if ( kk > end_x) then
1232 kk = end_x
1233 end if
1234 loop_2 : do l = 1,4
1235 stl(k,l) = 0.
1236 ll = j + l - 2
1237 if ( ll < start_y ) then
1238 ll = start_y
1239 else if ( ll > end_y) then
1240 ll = end_y
1241 end if
1242 stl(k,l) = array(kk,ll,izz)
1243 n = n + 1
1244 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1245 if (mask_relational == '>' .and. mask_array(kk,ll) > maskval) then
1246 is_masked = .true.
1247 else if (mask_relational == '<' .and. mask_array(kk,ll) < maskval) then
1248 is_masked = .true.
1249 else if (mask_relational == ' ' .and. mask_array(kk,ll) == maskval) then
1250 is_masked = .true.
1251 end if
1252 else if (present(mask_array) .and. present(maskval)) then
1253 if (mask_array(kk,ll) == maskval) is_masked = .true.
1254 end if
1255 if ( stl(k,l) == 0. .and. msgval /= 0.) then
1256 stl(k,l) = 1.E-20
1257 end if
1258 end do loop_2
1259 end do loop_1
1261 ! If we have a missing value, try the next interpolation method in the sequence
1262 if (present(mask_array) .and. present(maskval)) then
1263 do k=1,4
1264 do l=1,4
1265 if (stl(k,l) == msgval .or. is_masked) then
1266 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1267 start_z, end_z, msgval, interp_list, interp_opts, idx, &
1268 mask_relational, maskval, mask_array)
1269 return
1270 end if
1271 end do
1272 end do
1273 else
1274 do k=1,4
1275 do l=1,4
1276 if (stl(k,l) == msgval) then
1277 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, &
1278 start_z, end_z, msgval, interp_list, interp_opts, idx)
1279 return
1280 end if
1281 end do
1282 end do
1283 end if
1285 a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
1286 b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
1287 c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
1288 d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
1289 sixteen_pt = oned(y,a,b,c,d)
1291 if (n /= 16) then
1292 e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
1293 f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
1294 g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
1295 h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
1296 sixteen_pt = (sixteen_pt+oned(x,e,f,g,h)) * 0.5
1297 end if
1299 if (sixteen_pt == 1.E-20) sixteen_pt = 0.
1301 else
1302 if (present(mask_array) .and. present(maskval) .and. present(mask_relational)) then
1303 if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1304 mask_relational == '<' .and. mask_array(i,j) >= maskval .and. array(i,j,izz) /= msgval) then
1305 sixteen_pt = array(i,j,izz)
1306 else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1307 mask_relational == '>' .and. mask_array(i,j) <= maskval .and. array(i,j,izz) /= msgval) then
1308 sixteen_pt = array(i,j,izz)
1309 else if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1310 mask_relational == ' ' .and. mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
1311 sixteen_pt = array(i,j,izz)
1312 else
1313 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1314 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1315 end if
1316 else if (present(mask_array) .and. present(maskval)) then
1317 if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. &
1318 mask_array(i,j) /= maskval .and. array(i,j,izz) /= msgval) then
1319 sixteen_pt = array(i,j,izz)
1320 else
1321 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1322 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1323 end if
1324 else
1325 if (i >= start_x .and. i <= end_x .and. j >= start_y .and. j <= end_y .and. array(i,j,izz) /= msgval) then
1326 sixteen_pt = array(i,j,izz)
1327 else
1328 sixteen_pt = interp_sequence(xx, yy, izz, array, start_x, end_x, start_y, end_y, start_z, end_z, &
1329 msgval, interp_list, interp_opts, idx, mask_relational, maskval, mask_array)
1330 end if
1331 end if
1332 end if
1334 end function sixteen_pt
1337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1338 ! Name: oned
1339 !
1340 ! Purpose: 1-dimensional overlapping parabolic interpolation
1341 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342 function oned(x,a,b,c,d)
1344 implicit none
1346 ! Arguments
1347 real, intent(in) :: x,a,b,c,d
1349 ! Return value
1350 real :: oned
1352 oned = 0.
1354 if ( x == 0. ) then
1355 oned = b
1356 else if ( x == 1. ) then
1357 oned = c
1358 end if
1360 if (b*c /= 0.) then
1361 if ( a*d == 0. ) then
1362 if ( ( a == 0 ) .and. ( d == 0 ) ) then
1363 oned = b*(1.0-x)+c*x
1364 else if ( a /= 0. ) then
1365 oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))
1366 else if ( d /= 0. ) then
1367 oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c))
1368 end if
1369 else
1370 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)))
1371 end if
1372 end if
1374 end function oned
1376 end module interp_module