source: LMDZ5/branches/testing/libf/misc/regr1_conserv_m.F90 @ 2791

Last change on this file since 2791 was 2471, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2434:2457 into testing branch

File size: 10.6 KB
Line 
1module regr1_conserv_m
2
3  ! Author: Lionel GUEZ
4
5  use assert_eq_m, only: assert_eq
6  use assert_m, only: assert
7  use interpolation, only: locate
8
9  implicit none
10
11  interface regr1_conserv
12
13     ! This generic procedure regrids a piecewise linear function (not
14     ! necessarily continuous) by averaging it. This is a conservative
15     ! regridding. The regridding operation is done on the first
16     ! dimension of the input array. Input are positions of cell
17     ! edges. The target grid should be included in the source grid:
18     ! no extrapolation is allowed.
19
20     ! The only difference between the procedures is the rank of the
21     ! first argument.
22
23     ! real, intent(in), rank >= 1:: vs ! (ns, ...)
24     ! averages of cells of the source grid
25     ! vs(is, ...) for [xs(is), xs(is + 1)]
26
27     ! real, intent(in):: xs(:) ! (ns + 1)
28     ! edges of cells of the source grid, in strictly ascending order
29
30     ! real, intent(in):: xt(:) ! (nt + 1)
31     ! edges of cells of the target grid, in strictly ascending order
32
33     ! real, intent(in), optional, rank >= 1:: slope ! (ns, ...)
34     ! same rank as vs
35     ! slopes inside cells of the source grid
36     ! slope(is, ...) for [xs(is), xs(is + 1)]
37     ! If not present, slopes are taken equal to 0. The regridding is
38     ! then first order.
39
40     ! real, intent(out), rank >= 1:: vt(nt, ...)
41     ! has the same rank as vs and slope
42     ! averages of cells of the target grid
43     ! vt(it, ...) for  [xt(it), xt(it + 1)]
44
45     ! ns and nt must be >= 1.
46
47     ! See notes for explanations on the algorithm and justification
48     ! of algorithmic choices.
49
50     module procedure regr11_conserv, regr12_conserv, regr13_conserv, &
51          regr14_conserv
52  end interface regr1_conserv
53
54  private
55  public regr1_conserv
56
57contains
58
59  subroutine regr11_conserv(vs, xs, xt, vt, slope)
60
61    ! vs and slope have rank 1.
62
63    real, intent(in):: vs(:)
64    real, intent(in):: xs(:)
65    real, intent(in):: xt(:)
66    real, intent(out):: vt(:)
67    real, intent(in), optional:: slope(:)
68
69    ! Local:
70    integer is, it, ns, nt
71    logical slope_present
72
73    !---------------------------------------------
74
75    ns = assert_eq(size(vs), size(xs) - 1, "regr11_conserv ns")
76    nt = assert_eq(size(xt) - 1, size(vt), "regr11_conserv nt")
77
78    ! Quick check on sort order:
79    call assert(xs(1) < xs(2), "regr11_conserv xs bad order")
80    call assert(xt(1) < xt(2), "regr11_conserv xt bad order")
81
82    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
83         "regr11_conserv extrapolation")
84    slope_present = present(slope)
85
86    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
87    do it = 1, nt
88       ! 1 <= is <= ns
89       ! xs(is) <= xt(it) < xs(is + 1)
90       if (xt(it + 1) <= xs(is + 1)) then
91          vt(it) = mean_lin(xt(it), xt(it + 1))
92       else
93          vt(it) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
94          is = is + 1
95          do while (xs(is + 1) < xt(it + 1))
96             ! 1 <= is <= ns - 1
97             vt(it) = vt(it) + (xs(is + 1) - xs(is)) * vs(is)
98             is = is + 1
99          end do
100          ! 1 <= is <= ns
101          vt(it) = (vt(it) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) &
102               - xs(is))) / (xt(it + 1) - xt(it))
103       end if
104
105       if (xs(is + 1) == xt(it + 1)) is = is + 1
106       ! 1 <= is <= ns .or. it == nt
107    end do
108
109  contains
110
111    real function mean_lin(a, b)
112
113      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
114
115      real, intent(in):: a, b
116
117      !---------------------------------------------
118
119      if (slope_present) then
120         mean_lin = slope(is) / 2. * (a + b - xs(is) - xs(is + 1)) + vs(is)
121      else
122         mean_lin = vs(is)
123      end if
124
125    end function mean_lin
126
127  end subroutine regr11_conserv
128
129  !********************************************
130
131  subroutine regr12_conserv(vs, xs, xt, vt, slope)
132
133    ! vs and slope have rank 2.
134
135    real, intent(in):: vs(:, :)
136    real, intent(in):: xs(:)
137    real, intent(in):: xt(:)
138    real, intent(out):: vt(:, :)
139    real, intent(in), optional:: slope(:, :)
140
141    ! Local:
142    integer is, it, ns, nt, n2
143    logical slope_present
144
145    !---------------------------------------------
146
147    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_conserv ns")
148    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr12_conserv nt")
149    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr12_conserv n2")
150
151    ! Quick check on sort order:
152    call assert(xs(1) < xs(2), "regr12_conserv xs bad order")
153    call assert(xt(1) < xt(2), "regr12_conserv xt bad order")
154
155    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
156         "regr12_conserv extrapolation")
157    slope_present = present(slope)
158
159    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
160    do it = 1, nt
161       ! 1 <= is <= ns
162       ! xs(is) <= xt(it) < xs(is + 1)
163       if (xt(it + 1) <= xs(is + 1)) then
164          vt(it, :) = mean_lin(xt(it), xt(it + 1))
165       else
166          vt(it, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
167          is = is + 1
168          do while (xs(is + 1) < xt(it + 1))
169             ! 1 <= is <= ns - 1
170             vt(it, :) = vt(it, :) + (xs(is + 1) - xs(is)) * vs(is, :)
171             is = is + 1
172          end do
173          ! 1 <= is <= ns
174          vt(it, :) = (vt(it, :) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) &
175               - xs(is))) / (xt(it + 1) - xt(it))
176       end if
177
178       if (xs(is + 1) == xt(it + 1)) is = is + 1
179       ! 1 <= is <= ns .or. it == nt
180    end do
181
182  contains
183
184    function mean_lin(a, b)
185
186      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
187
188      real, intent(in):: a, b
189      real mean_lin(n2)
190
191      !---------------------------------------------
192
193      if (slope_present) then
194         mean_lin = slope(is, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
195              + vs(is, :)
196      else
197         mean_lin = vs(is, :)
198      end if
199
200    end function mean_lin
201
202  end subroutine regr12_conserv
203
204  !********************************************
205
206  subroutine regr13_conserv(vs, xs, xt, vt, slope)
207
208    ! vs and slope have rank 3.
209
210    real, intent(in):: vs(:, :, :)
211    real, intent(in):: xs(:)
212    real, intent(in):: xt(:)
213    real, intent(out):: vt(:, :, :)
214    real, intent(in), optional:: slope(:, :, :)
215
216    ! Local:
217    integer is, it, ns, nt, n2, n3
218    logical slope_present
219
220    !---------------------------------------------
221
222    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_conserv ns")
223    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr13_conserv nt")
224    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr13_conserv n2")
225    n3 = assert_eq(size(vs, 3), size(vt, 3), "regr13_conserv n3")
226
227    ! Quick check on sort order:
228    call assert(xs(1) < xs(2), "regr13_conserv xs bad order")
229    call assert(xt(1) < xt(2), "regr13_conserv xt bad order")
230
231    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
232         "regr13_conserv extrapolation")
233    slope_present = present(slope)
234
235    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
236    do it = 1, nt
237       ! 1 <= is <= ns
238       ! xs(is) <= xt(it) < xs(is + 1)
239       if (xt(it + 1) <= xs(is + 1)) then
240          vt(it, :, :) = mean_lin(xt(it), xt(it + 1))
241       else
242          vt(it, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
243          is = is + 1
244          do while (xs(is + 1) < xt(it + 1))
245             ! 1 <= is <= ns - 1
246             vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - xs(is)) * vs(is, :, :)
247             is = is + 1
248          end do
249          ! 1 <= is <= ns
250          vt(it, :, :) = (vt(it, :, :) + mean_lin(xs(is), xt(it + 1)) &
251               * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it))
252       end if
253
254       if (xs(is + 1) == xt(it + 1)) is = is + 1
255       ! 1 <= is <= ns .or. it == nt
256    end do
257
258  contains
259
260    function mean_lin(a, b)
261
262      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
263
264      real, intent(in):: a, b
265      real mean_lin(n2, n3)
266
267      !---------------------------------------------
268
269      if (slope_present) then
270         mean_lin = slope(is, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
271              + vs(is, :, :)
272      else
273         mean_lin = vs(is, :, :)
274      end if
275
276    end function mean_lin
277
278  end subroutine regr13_conserv
279
280  !********************************************
281
282  subroutine regr14_conserv(vs, xs, xt, vt, slope)
283
284    ! vs and slope have rank 4.
285
286    real, intent(in):: vs(:, :, :, :)
287    real, intent(in):: xs(:)
288    real, intent(in):: xt(:)
289    real, intent(out):: vt(:, :, :, :)
290    real, intent(in), optional:: slope(:, :, :, :)
291
292    ! Local:
293    integer is, it, ns, nt, n2, n3, n4
294    logical slope_present
295
296    !---------------------------------------------
297
298    ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_conserv ns")
299    nt = assert_eq(size(xt) - 1, size(vt, 1), "regr14_conserv nt")
300    n2 = assert_eq(size(vs, 2), size(vt, 2), "regr14_conserv n2")
301    n3 = assert_eq(size(vs, 3), size(vt, 3), "regr14_conserv n3")
302    n4 = assert_eq(size(vs, 4), size(vt, 4), "regr14_conserv n4")
303
304    ! Quick check on sort order:
305    call assert(xs(1) < xs(2), "regr14_conserv xs bad order")
306    call assert(xt(1) < xt(2), "regr14_conserv xt bad order")
307
308    call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), &
309         "regr14_conserv extrapolation")
310    slope_present = present(slope)
311
312    is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation
313    do it = 1, nt
314       ! 1 <= is <= ns
315       ! xs(is) <= xt(it) < xs(is + 1)
316       if (xt(it + 1) <= xs(is + 1)) then
317          vt(it, :, :, :) = mean_lin(xt(it), xt(it + 1))
318       else
319          vt(it, :, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it))
320          is = is + 1
321          do while (xs(is + 1) < xt(it + 1))
322             ! 1 <= is <= ns - 1
323             vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - xs(is)) &
324                  * vs(is, :, :, :)
325             is = is + 1
326          end do
327          ! 1 <= is <= ns
328          vt(it, :, :, :) = (vt(it, :, :, :) + mean_lin(xs(is), xt(it + 1)) &
329               * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it))
330       end if
331
332       if (xs(is + 1) == xt(it + 1)) is = is + 1
333       ! 1 <= is <= ns .or. it == nt
334    end do
335
336  contains
337
338    function mean_lin(a, b)
339
340      ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)]
341
342      real, intent(in):: a, b
343      real mean_lin(n2, n3, n4)
344
345      !---------------------------------------------
346
347      if (slope_present) then
348         mean_lin = slope(is, :, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) &
349              + vs(is, :, :, :)
350      else
351         mean_lin = vs(is, :, :, :)
352      end if
353
354    end function mean_lin
355
356  end subroutine regr14_conserv
357
358end module regr1_conserv_m
Note: See TracBrowser for help on using the repository browser.