trunk/libf/regr1_step_av.f90

module regr1_step_av_m

  ! Author: Lionel GUEZ

  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, &
          regr14_step_av
  end interface

  private
  public regr1_step_av

contains

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

    ! "vs" has rank 1.

    use nr_util, only: assert_eq, assert
    use numer_rec, 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 nr_util, only: assert_eq, assert
    use numer_rec, 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 nr_util, only: assert_eq, assert
    use numer_rec, 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

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

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

    ! "vs" has rank 4.

    use nr_util, only: assert_eq, assert
    use numer_rec, 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), size(vs, 4))
    ! (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, "regr14_step_av ns")
    nt = size(xt) - 1

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

    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
         "regr14_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 regr14_step_av

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