trunk/libf/regr3_lint.f90

module regr3_lint_m

  implicit none

  interface regr3_lint
     ! Each procedure regrids by linear interpolation.
     ! The regridding operation is done on the third dimension of the
     ! input array.
     ! The difference betwwen the procedures is the rank of the first argument.
     module procedure regr33_lint
  end interface

  private
  public regr3_lint

contains

  function regr33_lint(vs, xs, xt) result(vt)

    ! "vs" has rank 3.

    use numer_rec, only: assert_eq, hunt

    real, intent(in):: vs(:, :, :)
    ! (values of the function at source points "xs")

    real, intent(in):: xs(:)
    ! (abscissas of points in source grid, in strictly monotonic order)

    real, intent(in):: xt(:)
    ! (abscissas of points in target grid)

    real vt(size(vs, 1), size(vs, 2), size(xt))
    ! (values of the function on the target grid)

    ! Variables local to the procedure:
    integer is, it, ns
    integer is_b ! "is" bound between 1 and "ns - 1"

    !--------------------------------------

    ns = assert_eq(size(vs, 3), size(xs), "regr33_lint ns")

    is = -1 ! go immediately to bisection on first call to "hunt"

    do it = 1, size(xt)
       call hunt(xs, xt(it), is)
       is_b = min(max(is, 1), ns - 1)
       vt(:, :, it) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b) &
            + (xt(it) - xs(is_b)) * vs(:, :, is_b+1)) / (xs(is_b+1) - xs(is_b))
    end do

  end function regr33_lint

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