Changeset 1263 for LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90
- Timestamp:
- Nov 17, 2009, 2:00:14 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4-dev/libf/bibio/regr1_step_av_m.F90
r1157 r1263 17 17 ! The difference between the procedures is the rank of the first argument. 18 18 19 module procedure regr11_step_av, regr12_step_av, regr13_step_av 19 module procedure regr11_step_av, regr12_step_av, regr13_step_av, & 20 regr14_step_av 20 21 end interface 21 22 … … 203 204 end function regr13_step_av 204 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 205 268 end module regr1_step_av_m
Note: See TracChangeset
for help on using the changeset viewer.