| [1157] | 1 | ! $Id$ |
|---|
| 2 | module regr3_lint_m |
|---|
| 3 | |
|---|
| 4 | ! Author: Lionel GUEZ |
|---|
| 5 | |
|---|
| 6 | implicit none |
|---|
| 7 | |
|---|
| 8 | interface regr3_lint |
|---|
| 9 | ! Each procedure regrids by linear interpolation. |
|---|
| 10 | ! The regridding operation is done on the third dimension of the |
|---|
| 11 | ! input array. |
|---|
| 12 | ! The difference betwwen the procedures is the rank of the first argument. |
|---|
| 13 | module procedure regr33_lint |
|---|
| 14 | end interface |
|---|
| 15 | |
|---|
| 16 | private |
|---|
| 17 | public regr3_lint |
|---|
| 18 | |
|---|
| 19 | contains |
|---|
| 20 | |
|---|
| 21 | function regr33_lint(vs, xs, xt) result(vt) |
|---|
| 22 | |
|---|
| 23 | ! "vs" has rank 3. |
|---|
| 24 | |
|---|
| 25 | use assert_eq_m, only: assert_eq |
|---|
| 26 | use interpolation, only: hunt |
|---|
| 27 | |
|---|
| 28 | real, intent(in):: vs(:, :, :) |
|---|
| 29 | ! (values of the function at source points "xs") |
|---|
| 30 | |
|---|
| 31 | real, intent(in):: xs(:) |
|---|
| 32 | ! (abscissas of points in source grid, in strictly monotonic order) |
|---|
| 33 | |
|---|
| 34 | real, intent(in):: xt(:) |
|---|
| 35 | ! (abscissas of points in target grid) |
|---|
| 36 | |
|---|
| 37 | real vt(size(vs, 1), size(vs, 2), size(xt)) |
|---|
| 38 | ! (values of the function on the target grid) |
|---|
| 39 | |
|---|
| 40 | ! Variables local to the procedure: |
|---|
| 41 | integer is, it, ns |
|---|
| 42 | integer is_b ! "is" bound between 1 and "ns - 1" |
|---|
| 43 | |
|---|
| 44 | !-------------------------------------- |
|---|
| 45 | |
|---|
| 46 | ns = assert_eq(size(vs, 3), size(xs), "regr33_lint ns") |
|---|
| 47 | |
|---|
| 48 | is = -1 ! go immediately to bisection on first call to "hunt" |
|---|
| 49 | |
|---|
| 50 | do it = 1, size(xt) |
|---|
| 51 | call hunt(xs, xt(it), is) |
|---|
| 52 | is_b = min(max(is, 1), ns - 1) |
|---|
| 53 | vt(:, :, it) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b) & |
|---|
| 54 | + (xt(it) - xs(is_b)) * vs(:, :, is_b+1)) / (xs(is_b+1) - xs(is_b)) |
|---|
| 55 | end do |
|---|
| 56 | |
|---|
| 57 | end function regr33_lint |
|---|
| 58 | |
|---|
| 59 | end module regr3_lint_m |
|---|