1 |
module regr1_lint_m |
2 |
|
3 |
! Author: Lionel GUEZ |
4 |
|
5 |
implicit none |
6 |
|
7 |
interface regr1_lint |
8 |
! Each procedure regrids by linear interpolation. |
9 |
! The regridding operation is done on the first dimension of the |
10 |
! input array. |
11 |
! The difference betwwen the procedures is the rank of the first argument. |
12 |
module procedure regr11_lint, regr12_lint |
13 |
end interface |
14 |
|
15 |
private |
16 |
public regr1_lint |
17 |
|
18 |
contains |
19 |
|
20 |
function regr11_lint(vs, xs, xt) result(vt) |
21 |
|
22 |
! "vs" has rank 1. |
23 |
|
24 |
use nr_util, only: assert_eq |
25 |
use numer_rec, only: hunt !!, polint |
26 |
|
27 |
real, intent(in):: vs(:) |
28 |
! (values of the function at source points "xs") |
29 |
|
30 |
real, intent(in):: xs(:) |
31 |
! (abscissas of points in source grid, in strictly monotonic order) |
32 |
|
33 |
real, intent(in):: xt(:) |
34 |
! (abscissas of points in target grid) |
35 |
|
36 |
real vt(size(xt)) ! values of the function on the target grid |
37 |
|
38 |
! Variables local to the procedure: |
39 |
integer is, it, ns |
40 |
integer is_b ! "is" bound between 1 and "ns - 1" |
41 |
|
42 |
!-------------------------------------- |
43 |
|
44 |
ns = assert_eq(size(vs), size(xs), "regr11_lint ns") |
45 |
|
46 |
is = -1 ! go immediately to bisection on first call to "hunt" |
47 |
|
48 |
do it = 1, size(xt) |
49 |
call hunt(xs, xt(it), is) |
50 |
is_b = min(max(is, 1), ns - 1) |
51 |
!! call polint(xs(is_b:is_b+1), vs(is_b:is_b+1), xt(it), vt(it)) |
52 |
vt(it) = ((xs(is_b+1) - xt(it)) * vs(is_b) & |
53 |
+ (xt(it) - xs(is_b)) * vs(is_b+1)) / (xs(is_b+1) - xs(is_b)) |
54 |
end do |
55 |
|
56 |
end function regr11_lint |
57 |
|
58 |
!********************************************************* |
59 |
|
60 |
function regr12_lint(vs, xs, xt) result(vt) |
61 |
|
62 |
! "vs" has rank 2. |
63 |
|
64 |
use nr_util, only: assert_eq |
65 |
use numer_rec, only: hunt |
66 |
|
67 |
real, intent(in):: vs(:, :) |
68 |
! (values of the function at source points "xs") |
69 |
|
70 |
real, intent(in):: xs(:) |
71 |
! (abscissas of points in source grid, in strictly monotonic order) |
72 |
|
73 |
real, intent(in):: xt(:) |
74 |
! (abscissas of points in target grid) |
75 |
|
76 |
real vt(size(xt), size(vs, 2)) ! values of the function on the target grid |
77 |
|
78 |
! Variables local to the procedure: |
79 |
integer is, it, ns |
80 |
integer is_b ! "is" bound between 1 and "ns - 1" |
81 |
|
82 |
!-------------------------------------- |
83 |
|
84 |
ns = assert_eq(size(vs, 1), size(xs), "regr12_lint ns") |
85 |
|
86 |
is = -1 ! go immediately to bisection on first call to "hunt" |
87 |
|
88 |
do it = 1, size(xt) |
89 |
call hunt(xs, xt(it), is) |
90 |
is_b = min(max(is, 1), ns - 1) |
91 |
vt(it, :) = ((xs(is_b+1) - xt(it)) * vs(is_b, :) & |
92 |
+ (xt(it) - xs(is_b)) * vs(is_b+1, :)) / (xs(is_b+1) - xs(is_b)) |
93 |
end do |
94 |
|
95 |
end function regr12_lint |
96 |
|
97 |
end module regr1_lint_m |