| 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, regr34_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 | !********************************************************* |
|---|
| 60 | |
|---|
| 61 | function regr34_lint(vs, xs, xt) result(vt) |
|---|
| 62 | |
|---|
| 63 | ! "vs" has rank 4. |
|---|
| 64 | |
|---|
| 65 | use assert_eq_m, only: assert_eq |
|---|
| 66 | use interpolation, only: hunt |
|---|
| 67 | |
|---|
| 68 | real, intent(in):: vs(:, :, :, :) |
|---|
| 69 | ! (values of the function at source points "xs") |
|---|
| 70 | |
|---|
| 71 | real, intent(in):: xs(:) |
|---|
| 72 | ! (abscissas of points in source grid, in strictly monotonic order) |
|---|
| 73 | |
|---|
| 74 | real, intent(in):: xt(:) |
|---|
| 75 | ! (abscissas of points in target grid) |
|---|
| 76 | |
|---|
| 77 | real vt(size(vs, 1), size(vs, 2), size(xt), size(vs, 4)) |
|---|
| 78 | ! (values of the function on the target grid) |
|---|
| 79 | |
|---|
| 80 | ! Variables local to the procedure: |
|---|
| 81 | integer is, it, ns |
|---|
| 82 | integer is_b ! "is" bound between 1 and "ns - 1" |
|---|
| 83 | |
|---|
| 84 | !-------------------------------------- |
|---|
| 85 | |
|---|
| 86 | ns = assert_eq(size(vs, 3), size(xs), "regr34_lint ns") |
|---|
| 87 | |
|---|
| 88 | is = -1 ! go immediately to bisection on first call to "hunt" |
|---|
| 89 | |
|---|
| 90 | do it = 1, size(xt) |
|---|
| 91 | call hunt(xs, xt(it), is) |
|---|
| 92 | is_b = min(max(is, 1), ns - 1) |
|---|
| 93 | vt(:, :, it, :) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b, :) & |
|---|
| 94 | + (xt(it) - xs(is_b)) * vs(:, :, is_b+1, :)) & |
|---|
| 95 | / (xs(is_b+1) - xs(is_b)) |
|---|
| 96 | end do |
|---|
| 97 | |
|---|
| 98 | end function regr34_lint |
|---|
| 99 | |
|---|
| 100 | end module regr3_lint_m |
|---|