source: LMDZ4/branches/LMDZ4V5.0-dev/libf/bibio/regr1_step_av_m.F90 @ 1356

Last change on this file since 1356 was 1279, checked in by Laurent Fairhead, 15 years ago

Merged LMDZ4-dev branch changes r1241:1278 into the trunk
Running trunk and LMDZ4-dev in LMDZOR configuration on local
machine (sequential) and SX8 (4-proc) yields identical results
(restart and restartphy are identical binarily)
Log history from r1241 to r1278 is available by switching to
source:LMDZ4/branches/LMDZ4-dev-20091210

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