[1157] | 1 | ! $Id$ |
---|
| 2 | module 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 | 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 | use assert_eq_m, only: assert_eq |
---|
| 32 | use assert_m, only: assert |
---|
| 33 | use interpolation, only: locate |
---|
| 34 | |
---|
| 35 | real, intent(in):: vs(:) ! values of steps on the source grid |
---|
| 36 | ! (Step "is" is between "xs(is)" and "xs(is + 1)".) |
---|
| 37 | |
---|
| 38 | real, intent(in):: xs(:) |
---|
| 39 | ! (edges of of steps on the source grid, in strictly increasing order) |
---|
| 40 | |
---|
| 41 | real, intent(in):: xt(:) |
---|
| 42 | ! (edges of cells of the target grid, in strictly increasing order) |
---|
| 43 | |
---|
| 44 | real vt(size(xt) - 1) ! average values on the target grid |
---|
| 45 | ! (Cell "it" is between "xt(it)" and "xt(it + 1)".) |
---|
| 46 | |
---|
| 47 | ! Variables local to the procedure: |
---|
| 48 | integer is, it, ns, nt |
---|
| 49 | real left_edge |
---|
| 50 | |
---|
| 51 | !--------------------------------------------- |
---|
| 52 | |
---|
| 53 | ns = assert_eq(size(vs), size(xs) - 1, "regr11_step_av ns") |
---|
| 54 | nt = size(xt) - 1 |
---|
| 55 | ! Quick check on sort order: |
---|
| 56 | call assert(xs(1) < xs(2), "regr11_step_av xs bad order") |
---|
| 57 | call assert(xt(1) < xt(2), "regr11_step_av xt bad order") |
---|
| 58 | |
---|
| 59 | call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & |
---|
| 60 | "regr11_step_av extrapolation") |
---|
| 61 | |
---|
| 62 | is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation |
---|
| 63 | do it = 1, nt |
---|
| 64 | ! 1 <= is <= ns |
---|
| 65 | ! xs(is) <= xt(it) < xs(is + 1) |
---|
| 66 | ! Compute "vt(it)": |
---|
| 67 | left_edge = xt(it) |
---|
| 68 | vt(it) = 0. |
---|
| 69 | do while (xs(is + 1) < xt(it + 1)) |
---|
| 70 | ! 1 <= is <= ns - 1 |
---|
| 71 | vt(it) = vt(it) + (xs(is + 1) - left_edge) * vs(is) |
---|
| 72 | is = is + 1 |
---|
| 73 | left_edge = xs(is) |
---|
| 74 | end do |
---|
| 75 | ! 1 <= is <= ns |
---|
| 76 | vt(it) = (vt(it) + (xt(it + 1) - left_edge) * vs(is)) & |
---|
| 77 | / (xt(it + 1) - xt(it)) |
---|
| 78 | if (xs(is + 1) == xt(it + 1)) is = is + 1 |
---|
| 79 | ! 1 <= is <= ns .or. it == nt |
---|
| 80 | end do |
---|
| 81 | |
---|
| 82 | end function regr11_step_av |
---|
| 83 | |
---|
| 84 | !******************************************** |
---|
| 85 | |
---|
| 86 | function regr12_step_av(vs, xs, xt) result(vt) |
---|
| 87 | |
---|
| 88 | ! "vs" has rank 2. |
---|
| 89 | |
---|
| 90 | use assert_eq_m, only: assert_eq |
---|
| 91 | use assert_m, only: assert |
---|
| 92 | use interpolation, only: locate |
---|
| 93 | |
---|
| 94 | real, intent(in):: vs(:, :) ! values of steps on the source grid |
---|
| 95 | ! (Step "is" is between "xs(is)" and "xs(is + 1)".) |
---|
| 96 | |
---|
| 97 | real, intent(in):: xs(:) |
---|
| 98 | ! (edges of steps on the source grid, in strictly increasing order) |
---|
| 99 | |
---|
| 100 | real, intent(in):: xt(:) |
---|
| 101 | ! (edges of cells of the target grid, in strictly increasing order) |
---|
| 102 | |
---|
| 103 | real vt(size(xt) - 1, size(vs, 2)) ! average values on the target grid |
---|
| 104 | ! (Cell "it" is between "xt(it)" and "xt(it + 1)".) |
---|
| 105 | |
---|
| 106 | ! Variables local to the procedure: |
---|
| 107 | integer is, it, ns, nt |
---|
| 108 | real left_edge |
---|
| 109 | |
---|
| 110 | !--------------------------------------------- |
---|
| 111 | |
---|
| 112 | ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_step_av ns") |
---|
| 113 | nt = size(xt) - 1 |
---|
| 114 | |
---|
| 115 | ! Quick check on sort order: |
---|
| 116 | call assert(xs(1) < xs(2), "regr12_step_av xs bad order") |
---|
| 117 | call assert(xt(1) < xt(2), "regr12_step_av xt bad order") |
---|
| 118 | |
---|
| 119 | call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & |
---|
| 120 | "regr12_step_av extrapolation") |
---|
| 121 | |
---|
| 122 | is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation |
---|
| 123 | do it = 1, nt |
---|
| 124 | ! 1 <= is <= ns |
---|
| 125 | ! xs(is) <= xt(it) < xs(is + 1) |
---|
| 126 | ! Compute "vt(it, :)": |
---|
| 127 | left_edge = xt(it) |
---|
| 128 | vt(it, :) = 0. |
---|
| 129 | do while (xs(is + 1) < xt(it + 1)) |
---|
| 130 | ! 1 <= is <= ns - 1 |
---|
| 131 | vt(it, :) = vt(it, :) + (xs(is + 1) - left_edge) * vs(is, :) |
---|
| 132 | is = is + 1 |
---|
| 133 | left_edge = xs(is) |
---|
| 134 | end do |
---|
| 135 | ! 1 <= is <= ns |
---|
| 136 | vt(it, :) = (vt(it, :) + (xt(it + 1) - left_edge) * vs(is, :)) & |
---|
| 137 | / (xt(it + 1) - xt(it)) |
---|
| 138 | if (xs(is + 1) == xt(it + 1)) is = is + 1 |
---|
| 139 | ! 1 <= is <= ns .or. it == nt |
---|
| 140 | end do |
---|
| 141 | |
---|
| 142 | end function regr12_step_av |
---|
| 143 | |
---|
| 144 | !******************************************** |
---|
| 145 | |
---|
| 146 | function regr13_step_av(vs, xs, xt) result(vt) |
---|
| 147 | |
---|
| 148 | ! "vs" has rank 3. |
---|
| 149 | |
---|
| 150 | use assert_eq_m, only: assert_eq |
---|
| 151 | use assert_m, only: assert |
---|
| 152 | use interpolation, only: locate |
---|
| 153 | |
---|
| 154 | real, intent(in):: vs(:, :, :) ! values of steps on the source grid |
---|
| 155 | ! (Step "is" is between "xs(is)" and "xs(is + 1)".) |
---|
| 156 | |
---|
| 157 | real, intent(in):: xs(:) |
---|
| 158 | ! (edges of steps on the source grid, in strictly increasing order) |
---|
| 159 | |
---|
| 160 | real, intent(in):: xt(:) |
---|
| 161 | ! (edges of cells of the target grid, in strictly increasing order) |
---|
| 162 | |
---|
| 163 | real vt(size(xt) - 1, size(vs, 2), size(vs, 3)) |
---|
| 164 | ! (average values on the target grid) |
---|
| 165 | ! (Cell "it" is between "xt(it)" and "xt(it + 1)".) |
---|
| 166 | |
---|
| 167 | ! Variables local to the procedure: |
---|
| 168 | integer is, it, ns, nt |
---|
| 169 | real left_edge |
---|
| 170 | |
---|
| 171 | !--------------------------------------------- |
---|
| 172 | |
---|
| 173 | ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_step_av ns") |
---|
| 174 | nt = size(xt) - 1 |
---|
| 175 | |
---|
| 176 | ! Quick check on sort order: |
---|
| 177 | call assert(xs(1) < xs(2), "regr13_step_av xs bad order") |
---|
| 178 | call assert(xt(1) < xt(2), "regr13_step_av xt bad order") |
---|
| 179 | |
---|
| 180 | call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & |
---|
| 181 | "regr13_step_av extrapolation") |
---|
| 182 | |
---|
| 183 | is = locate(xs, xt(1)) ! 1 <= is <= ns, because we forbid extrapolation |
---|
| 184 | do it = 1, nt |
---|
| 185 | ! 1 <= is <= ns |
---|
| 186 | ! xs(is) <= xt(it) < xs(is + 1) |
---|
| 187 | ! Compute "vt(it, :, :)": |
---|
| 188 | left_edge = xt(it) |
---|
| 189 | vt(it, :, :) = 0. |
---|
| 190 | do while (xs(is + 1) < xt(it + 1)) |
---|
| 191 | ! 1 <= is <= ns - 1 |
---|
| 192 | vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - left_edge) * vs(is, :, :) |
---|
| 193 | is = is + 1 |
---|
| 194 | left_edge = xs(is) |
---|
| 195 | end do |
---|
| 196 | ! 1 <= is <= ns |
---|
| 197 | vt(it, :, :) = (vt(it, :, :) & |
---|
| 198 | + (xt(it + 1) - left_edge) * vs(is, :, :)) / (xt(it + 1) - xt(it)) |
---|
| 199 | if (xs(is + 1) == xt(it + 1)) is = is + 1 |
---|
| 200 | ! 1 <= is <= ns .or. it == nt |
---|
| 201 | end do |
---|
| 202 | |
---|
| 203 | end function regr13_step_av |
---|
| 204 | |
---|
| 205 | end module regr1_step_av_m |
---|