trunk/libf/regr1_step_av.f90

module regr1_step_av_m

  implicit none

  interface regr1_step_av

     ! Each procedure regrids a step function by averaging it.
     ! The regridding operation is done on the first dimension of the
     ! input array.
     ! Source grid contains edges of steps.
     ! Target grid contains positions of cell edges.
     ! The target grid should be included in the source grid: no
     ! extrapolation is allowed.
     ! The difference between the procedures is the rank of the first argument.

     module procedure regr11_step_av, regr12_step_av, regr13_step_av
  end interface

  private
  public regr1_step_av

contains

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

    ! "vs" has rank 1.

    use nrutil, only: assert_eq, assert
    use interpolation, only: locate

    real, intent(in):: vs(:) ! values of steps on the source grid
    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)

    real, intent(in):: xs(:)
    ! (edges of of steps on the source grid, in strictly increasing order)

    real, intent(in):: xt(:)
    ! (edges of cells of the target grid, in strictly increasing order)

    real vt(size(xt) - 1) ! average values on the target grid
    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)

    ! Variables local to the procedure:
    integer is, it, ns, nt
    real left_edge

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

    ns = assert_eq(size(vs), size(xs) - 1, "regr11_step_av ns")
    nt = size(xt) - 1
    ! Quick check on sort order:
    call assert(xs(1) < xs(2), "regr11_step_av xs bad order")
    call assert(xt(1) < xt(2), "regr11_step_av xt bad order")

    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
         "regr11_step_av extrapolation")

    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
    do it = 1, nt
       ! 1 <= is <= ns
       ! xs(is) <= xt(it) < xs(is + 1)
       ! Compute "vt(it)":
       left_edge = xt(it)
       vt(it) = 0.
       do while (xs(is + 1) < xt(it + 1))
          ! 1 <= is <= ns - 1
          vt(it) = vt(it) + (xs(is + 1) - left_edge) * vs(is)
          is = is + 1
          left_edge = xs(is)
       end do
       ! 1 <= is <= ns
       vt(it) = (vt(it) + (xt(it + 1) - left_edge) * vs(is)) &
            / (xt(it + 1) - xt(it))
       if (xs(is + 1) == xt(it + 1)) is = is + 1
       ! 1 <= is <= ns .or. it == nt
    end do

  end function regr11_step_av

  !********************************************

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

    ! "vs" has rank 2.

    use nrutil, only: assert_eq, assert
    use interpolation, only: locate

    real, intent(in):: vs(:, :) ! values of steps on the source grid
    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)

    real, intent(in):: xs(:)
    ! (edges of steps on the source grid, in strictly increasing order)

    real, intent(in):: xt(:)
    ! (edges of cells of the target grid, in strictly increasing order)

    real vt(size(xt) - 1, size(vs, 2)) ! average values on the target grid
    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)

    ! Variables local to the procedure:
    integer is, it, ns, nt
    real left_edge

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

    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_step_av ns")
    nt = size(xt) - 1

    ! Quick check on sort order:
    call assert(xs(1) < xs(2), "regr12_step_av xs bad order")
    call assert(xt(1) < xt(2), "regr12_step_av xt bad order")

    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
         "regr12_step_av extrapolation")

    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
    do it = 1, nt
       ! 1 <= is <= ns
       ! xs(is) <= xt(it) < xs(is + 1)
       ! Compute "vt(it, :)":
       left_edge = xt(it)
       vt(it, :) = 0.
       do while (xs(is + 1) < xt(it + 1))
          ! 1 <= is <= ns - 1
          vt(it, :) = vt(it, :) + (xs(is + 1) - left_edge) * vs(is, :)
          is = is + 1
          left_edge = xs(is)
       end do
       ! 1 <= is <= ns
       vt(it, :) = (vt(it, :) + (xt(it + 1) - left_edge) * vs(is, :)) &
            / (xt(it + 1) - xt(it))
       if (xs(is + 1) == xt(it + 1)) is = is + 1
       ! 1 <= is <= ns .or. it == nt
    end do

  end function regr12_step_av

  !********************************************

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

    ! "vs" has rank 3.

    use nrutil, only: assert_eq, assert
    use interpolation, only: locate

    real, intent(in):: vs(:, :, :) ! values of steps on the source grid
    ! (Step "is" is between "xs(is)" and "xs(is + 1)".)

    real, intent(in):: xs(:)
    ! (edges of steps on the source grid, in strictly increasing order)

    real, intent(in):: xt(:)
    ! (edges of cells of the target grid, in strictly increasing order)

    real vt(size(xt) - 1, size(vs, 2), size(vs, 3)) 
    ! (average values on the target grid)
    ! (Cell "it" is between "xt(it)" and "xt(it + 1)".)

    ! Variables local to the procedure:
    integer is, it, ns, nt
    real left_edge

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

    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_step_av ns")
    nt = size(xt) - 1

    ! Quick check on sort order:
    call assert(xs(1) < xs(2), "regr13_step_av xs bad order")
    call assert(xt(1) < xt(2), "regr13_step_av xt bad order")

    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
         "regr13_step_av extrapolation")

    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
    do it = 1, nt
       ! 1 <= is <= ns
       ! xs(is) <= xt(it) < xs(is + 1)
       ! Compute "vt(it, :, :)":
       left_edge = xt(it)
       vt(it, :, :) = 0.
       do while (xs(is + 1) < xt(it + 1))
          ! 1 <= is <= ns - 1
          vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - left_edge) * vs(is, :, :)
          is = is + 1
          left_edge = xs(is)
       end do
       ! 1 <= is <= ns
       vt(it, :, :) = (vt(it, :, :) &
            + (xt(it + 1) - left_edge) * vs(is, :, :)) / (xt(it + 1) - xt(it))
       if (xs(is + 1) == xt(it + 1)) is = is + 1
       ! 1 <= is <= ns .or. it == nt
    end do

  end function regr13_step_av

end module regr1_step_av_m
1	guez	9	module regr1_step_av_m
2
3			implicit none
4
5			interface regr1_step_av
6
7			! Each procedure regrids a step function by averaging it.
8			! The regridding operation is done on the first dimension of the
9			! input array.
10			! Source grid contains edges of steps.
11			! Target grid contains positions of cell edges.
12			! The target grid should be included in the source grid: no
13			! extrapolation is allowed.
14			! The difference between the procedures is the rank of the first argument.
15
16			module procedure regr11_step_av, regr12_step_av, regr13_step_av
17			end interface
18
19			private
20			public regr1_step_av
21
22			contains
23
24			function regr11_step_av(vs, xs, xt) result(vt)
25
26			! "vs" has rank 1.
27
28			use nrutil, only: assert_eq, assert
29			use interpolation, only: locate
30
31			real, intent(in):: vs(:) ! values of steps on the source grid
32			! (Step "is" is between "xs(is)" and "xs(is + 1)".)
33
34			real, intent(in):: xs(:)
35			! (edges of of steps on the source grid, in strictly increasing order)
36
37			real, intent(in):: xt(:)
38			! (edges of cells of the target grid, in strictly increasing order)
39
40			real vt(size(xt) - 1) ! average values on the target grid
41			! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
42
43			! Variables local to the procedure:
44			integer is, it, ns, nt
45			real left_edge
46
47			!---------------------------------------------
48
49			ns = assert_eq(size(vs), size(xs) - 1, "regr11_step_av ns")
50			nt = size(xt) - 1
51			! Quick check on sort order:
52			call assert(xs(1) < xs(2), "regr11_step_av xs bad order")
53			call assert(xt(1) < xt(2), "regr11_step_av xt bad order")
54
55			call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
56			"regr11_step_av extrapolation")
57
58			is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
59			do it = 1, nt
60			! 1 <= is <= ns
61			! xs(is) <= xt(it) < xs(is + 1)
62			! Compute "vt(it)":
63			left_edge = xt(it)
64			vt(it) = 0.
65			do while (xs(is + 1) < xt(it + 1))
66			! 1 <= is <= ns - 1
67			vt(it) = vt(it) + (xs(is + 1) - left_edge) * vs(is)
68			is = is + 1
69			left_edge = xs(is)
70			end do
71			! 1 <= is <= ns
72			vt(it) = (vt(it) + (xt(it + 1) - left_edge) * vs(is)) &
73			/ (xt(it + 1) - xt(it))
74			if (xs(is + 1) == xt(it + 1)) is = is + 1
75			! 1 <= is <= ns .or. it == nt
76			end do
77
78			end function regr11_step_av
79
80			!********************************************
81
82			function regr12_step_av(vs, xs, xt) result(vt)
83
84			! "vs" has rank 2.
85
86			use nrutil, only: assert_eq, assert
87			use interpolation, only: locate
88
89			real, intent(in):: vs(:, :) ! values of steps on the source grid
90			! (Step "is" is between "xs(is)" and "xs(is + 1)".)
91
92			real, intent(in):: xs(:)
93			! (edges of steps on the source grid, in strictly increasing order)
94
95			real, intent(in):: xt(:)
96			! (edges of cells of the target grid, in strictly increasing order)
97
98			real vt(size(xt) - 1, size(vs, 2)) ! average values on the target grid
99			! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
100
101			! Variables local to the procedure:
102			integer is, it, ns, nt
103			real left_edge
104
105			!---------------------------------------------
106
107			ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_step_av ns")
108			nt = size(xt) - 1
109
110			! Quick check on sort order:
111			call assert(xs(1) < xs(2), "regr12_step_av xs bad order")
112			call assert(xt(1) < xt(2), "regr12_step_av xt bad order")
113
114			call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
115			"regr12_step_av extrapolation")
116
117			is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
118			do it = 1, nt
119			! 1 <= is <= ns
120			! xs(is) <= xt(it) < xs(is + 1)
121			! Compute "vt(it, :)":
122			left_edge = xt(it)
123			vt(it, :) = 0.
124			do while (xs(is + 1) < xt(it + 1))
125			! 1 <= is <= ns - 1
126			vt(it, :) = vt(it, :) + (xs(is + 1) - left_edge) * vs(is, :)
127			is = is + 1
128			left_edge = xs(is)
129			end do
130			! 1 <= is <= ns
131			vt(it, :) = (vt(it, :) + (xt(it + 1) - left_edge) * vs(is, :)) &
132			/ (xt(it + 1) - xt(it))
133			if (xs(is + 1) == xt(it + 1)) is = is + 1
134			! 1 <= is <= ns .or. it == nt
135			end do
136
137			end function regr12_step_av
138
139			!********************************************
140
141			function regr13_step_av(vs, xs, xt) result(vt)
142
143			! "vs" has rank 3.
144
145			use nrutil, only: assert_eq, assert
146			use interpolation, only: locate
147
148			real, intent(in):: vs(:, :, :) ! values of steps on the source grid
149			! (Step "is" is between "xs(is)" and "xs(is + 1)".)
150
151			real, intent(in):: xs(:)
152			! (edges of steps on the source grid, in strictly increasing order)
153
154			real, intent(in):: xt(:)
155			! (edges of cells of the target grid, in strictly increasing order)
156
157			real vt(size(xt) - 1, size(vs, 2), size(vs, 3))
158			! (average values on the target grid)
159			! (Cell "it" is between "xt(it)" and "xt(it + 1)".)
160
161			! Variables local to the procedure:
162			integer is, it, ns, nt
163			real left_edge
164
165			!---------------------------------------------
166
167			ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_step_av ns")
168			nt = size(xt) - 1
169
170			! Quick check on sort order:
171			call assert(xs(1) < xs(2), "regr13_step_av xs bad order")
172			call assert(xt(1) < xt(2), "regr13_step_av xt bad order")
173
174			call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
175			"regr13_step_av extrapolation")
176
177			is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
178			do it = 1, nt
179			! 1 <= is <= ns
180			! xs(is) <= xt(it) < xs(is + 1)
181			! Compute "vt(it, :, :)":
182			left_edge = xt(it)
183			vt(it, :, :) = 0.
184			do while (xs(is + 1) < xt(it + 1))
185			! 1 <= is <= ns - 1
186			vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - left_edge) * vs(is, :, :)
187			is = is + 1
188			left_edge = xs(is)
189			end do
190			! 1 <= is <= ns
191			vt(it, :, :) = (vt(it, :, :) &
192			+ (xt(it + 1) - left_edge) * vs(is, :, :)) / (xt(it + 1) - xt(it))
193			if (xs(is + 1) == xt(it + 1)) is = is + 1
194			! 1 <= is <= ns .or. it == nt
195			end do
196
197			end function regr13_step_av
198
199			end module regr1_step_av_m