Index: trunk/LMDZ.COMMON/libf/misc/regr1_conserv_m.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/misc/regr1_conserv_m.F90 (revision 1592) +++ trunk/LMDZ.COMMON/libf/misc/regr1_conserv_m.F90 (revision 1592) @@ -0,0 +1,358 @@ +module regr1_conserv_m + + ! Author: Lionel GUEZ + + use assert_eq_m, only: assert_eq + use assert_m, only: assert + use interpolation, only: locate + + implicit none + + interface regr1_conserv + + ! This generic procedure regrids a piecewise linear function (not + ! necessarily continuous) by averaging it. This is a conservative + ! regridding. The regridding operation is done on the first + ! dimension of the input array. Input are positions of cell + ! edges. The target grid should be included in the source grid: + ! no extrapolation is allowed. + + ! The only difference between the procedures is the rank of the + ! first argument. + + ! real, intent(in), rank >= 1:: vs ! (ns, ...) + ! averages of cells of the source grid + ! vs(is, ...) for [xs(is), xs(is + 1)] + + ! real, intent(in):: xs(:) ! (ns + 1) + ! edges of cells of the source grid, in strictly ascending order + + ! real, intent(in):: xt(:) ! (nt + 1) + ! edges of cells of the target grid, in strictly ascending order + + ! real, intent(in), optional, rank >= 1:: slope ! (ns, ...) + ! same rank as vs + ! slopes inside cells of the source grid + ! slope(is, ...) for [xs(is), xs(is + 1)] + ! If not present, slopes are taken equal to 0. The regridding is + ! then first order. + + ! real, intent(out), rank >= 1:: vt(nt, ...) + ! has the same rank as vs and slope + ! averages of cells of the target grid + ! vt(it, ...) for [xt(it), xt(it + 1)] + + ! ns and nt must be >= 1. + + ! See notes for explanations on the algorithm and justification + ! of algorithmic choices. + + module procedure regr11_conserv, regr12_conserv, regr13_conserv, & + regr14_conserv + end interface regr1_conserv + + private + public regr1_conserv + +contains + + subroutine regr11_conserv(vs, xs, xt, vt, slope) + + ! vs and slope have rank 1. + + real, intent(in):: vs(:) + real, intent(in):: xs(:) + real, intent(in):: xt(:) + real, intent(out):: vt(:) + real, intent(in), optional:: slope(:) + + ! Local: + integer is, it, ns, nt + logical slope_present + + !--------------------------------------------- + + ns = assert_eq(size(vs), size(xs) - 1, "regr11_conserv ns") + nt = assert_eq(size(xt) - 1, size(vt), "regr11_conserv nt") + + ! Quick check on sort order: + call assert(xs(1) < xs(2), "regr11_conserv xs bad order") + call assert(xt(1) < xt(2), "regr11_conserv xt bad order") + + call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & + "regr11_conserv extrapolation") + slope_present = present(slope) + + 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) + if (xt(it + 1) <= xs(is + 1)) then + vt(it) = mean_lin(xt(it), xt(it + 1)) + else + vt(it) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it)) + is = is + 1 + do while (xs(is + 1) < xt(it + 1)) + ! 1 <= is <= ns - 1 + vt(it) = vt(it) + (xs(is + 1) - xs(is)) * vs(is) + is = is + 1 + end do + ! 1 <= is <= ns + vt(it) = (vt(it) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) & + - xs(is))) / (xt(it + 1) - xt(it)) + end if + + if (xs(is + 1) == xt(it + 1)) is = is + 1 + ! 1 <= is <= ns .or. it == nt + end do + + contains + + real function mean_lin(a, b) + + ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)] + + real, intent(in):: a, b + + !--------------------------------------------- + + if (slope_present) then + mean_lin = slope(is) / 2. * (a + b - xs(is) - xs(is + 1)) + vs(is) + else + mean_lin = vs(is) + end if + + end function mean_lin + + end subroutine regr11_conserv + + !******************************************** + + subroutine regr12_conserv(vs, xs, xt, vt, slope) + + ! vs and slope have rank 2. + + real, intent(in):: vs(:, :) + real, intent(in):: xs(:) + real, intent(in):: xt(:) + real, intent(out):: vt(:, :) + real, intent(in), optional:: slope(:, :) + + ! Local: + integer is, it, ns, nt, n2 + logical slope_present + + !--------------------------------------------- + + ns = assert_eq(size(vs, 1), size(xs) - 1, "regr12_conserv ns") + nt = assert_eq(size(xt) - 1, size(vt, 1), "regr12_conserv nt") + n2 = assert_eq(size(vs, 2), size(vt, 2), "regr12_conserv n2") + + ! Quick check on sort order: + call assert(xs(1) < xs(2), "regr12_conserv xs bad order") + call assert(xt(1) < xt(2), "regr12_conserv xt bad order") + + call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & + "regr12_conserv extrapolation") + slope_present = present(slope) + + 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) + if (xt(it + 1) <= xs(is + 1)) then + vt(it, :) = mean_lin(xt(it), xt(it + 1)) + else + vt(it, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it)) + is = is + 1 + do while (xs(is + 1) < xt(it + 1)) + ! 1 <= is <= ns - 1 + vt(it, :) = vt(it, :) + (xs(is + 1) - xs(is)) * vs(is, :) + is = is + 1 + end do + ! 1 <= is <= ns + vt(it, :) = (vt(it, :) + mean_lin(xs(is), xt(it + 1)) * (xt(it + 1) & + - xs(is))) / (xt(it + 1) - xt(it)) + end if + + if (xs(is + 1) == xt(it + 1)) is = is + 1 + ! 1 <= is <= ns .or. it == nt + end do + + contains + + function mean_lin(a, b) + + ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)] + + real, intent(in):: a, b + real mean_lin(n2) + + !--------------------------------------------- + + if (slope_present) then + mean_lin = slope(is, :) / 2. * (a + b - xs(is) - xs(is + 1)) & + + vs(is, :) + else + mean_lin = vs(is, :) + end if + + end function mean_lin + + end subroutine regr12_conserv + + !******************************************** + + subroutine regr13_conserv(vs, xs, xt, vt, slope) + + ! vs and slope have rank 3. + + real, intent(in):: vs(:, :, :) + real, intent(in):: xs(:) + real, intent(in):: xt(:) + real, intent(out):: vt(:, :, :) + real, intent(in), optional:: slope(:, :, :) + + ! Local: + integer is, it, ns, nt, n2, n3 + logical slope_present + + !--------------------------------------------- + + ns = assert_eq(size(vs, 1), size(xs) - 1, "regr13_conserv ns") + nt = assert_eq(size(xt) - 1, size(vt, 1), "regr13_conserv nt") + n2 = assert_eq(size(vs, 2), size(vt, 2), "regr13_conserv n2") + n3 = assert_eq(size(vs, 3), size(vt, 3), "regr13_conserv n3") + + ! Quick check on sort order: + call assert(xs(1) < xs(2), "regr13_conserv xs bad order") + call assert(xt(1) < xt(2), "regr13_conserv xt bad order") + + call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & + "regr13_conserv extrapolation") + slope_present = present(slope) + + 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) + if (xt(it + 1) <= xs(is + 1)) then + vt(it, :, :) = mean_lin(xt(it), xt(it + 1)) + else + vt(it, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it)) + is = is + 1 + do while (xs(is + 1) < xt(it + 1)) + ! 1 <= is <= ns - 1 + vt(it, :, :) = vt(it, :, :) + (xs(is + 1) - xs(is)) * vs(is, :, :) + is = is + 1 + end do + ! 1 <= is <= ns + vt(it, :, :) = (vt(it, :, :) + mean_lin(xs(is), xt(it + 1)) & + * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it)) + end if + + if (xs(is + 1) == xt(it + 1)) is = is + 1 + ! 1 <= is <= ns .or. it == nt + end do + + contains + + function mean_lin(a, b) + + ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)] + + real, intent(in):: a, b + real mean_lin(n2, n3) + + !--------------------------------------------- + + if (slope_present) then + mean_lin = slope(is, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) & + + vs(is, :, :) + else + mean_lin = vs(is, :, :) + end if + + end function mean_lin + + end subroutine regr13_conserv + + !******************************************** + + subroutine regr14_conserv(vs, xs, xt, vt, slope) + + ! vs and slope have rank 4. + + real, intent(in):: vs(:, :, :, :) + real, intent(in):: xs(:) + real, intent(in):: xt(:) + real, intent(out):: vt(:, :, :, :) + real, intent(in), optional:: slope(:, :, :, :) + + ! Local: + integer is, it, ns, nt, n2, n3, n4 + logical slope_present + + !--------------------------------------------- + + ns = assert_eq(size(vs, 1), size(xs) - 1, "regr14_conserv ns") + nt = assert_eq(size(xt) - 1, size(vt, 1), "regr14_conserv nt") + n2 = assert_eq(size(vs, 2), size(vt, 2), "regr14_conserv n2") + n3 = assert_eq(size(vs, 3), size(vt, 3), "regr14_conserv n3") + n4 = assert_eq(size(vs, 4), size(vt, 4), "regr14_conserv n4") + + ! Quick check on sort order: + call assert(xs(1) < xs(2), "regr14_conserv xs bad order") + call assert(xt(1) < xt(2), "regr14_conserv xt bad order") + + call assert(xs(1) <= xt(1) .and. xt(nt + 1) <= xs(ns + 1), & + "regr14_conserv extrapolation") + slope_present = present(slope) + + 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) + if (xt(it + 1) <= xs(is + 1)) then + vt(it, :, :, :) = mean_lin(xt(it), xt(it + 1)) + else + vt(it, :, :, :) = mean_lin(xt(it), xs(is + 1)) * (xs(is + 1) - xt(it)) + is = is + 1 + do while (xs(is + 1) < xt(it + 1)) + ! 1 <= is <= ns - 1 + vt(it, :, :, :) = vt(it, :, :, :) + (xs(is + 1) - xs(is)) & + * vs(is, :, :, :) + is = is + 1 + end do + ! 1 <= is <= ns + vt(it, :, :, :) = (vt(it, :, :, :) + mean_lin(xs(is), xt(it + 1)) & + * (xt(it + 1) - xs(is))) / (xt(it + 1) - xt(it)) + end if + + if (xs(is + 1) == xt(it + 1)) is = is + 1 + ! 1 <= is <= ns .or. it == nt + end do + + contains + + function mean_lin(a, b) + + ! mean in [a, b] of the linear function in [xs(is), xs(is + 1)] + + real, intent(in):: a, b + real mean_lin(n2, n3, n4) + + !--------------------------------------------- + + if (slope_present) then + mean_lin = slope(is, :, :, :) / 2. * (a + b - xs(is) - xs(is + 1)) & + + vs(is, :, :, :) + else + mean_lin = vs(is, :, :, :) + end if + + end function mean_lin + + end subroutine regr14_conserv + +end module regr1_conserv_m Index: trunk/LMDZ.COMMON/libf/misc/slopes_m.F90 =================================================================== --- trunk/LMDZ.COMMON/libf/misc/slopes_m.F90 (revision 1592) +++ trunk/LMDZ.COMMON/libf/misc/slopes_m.F90 (revision 1592) @@ -0,0 +1,226 @@ +module slopes_m + + ! Author: Lionel GUEZ + + implicit none + + interface slopes + ! This generic function computes second order slopes with Van + ! Leer slope-limiting, given cell averages. Reference: Dukowicz, + ! 1987, SIAM Journal on Scientific and Statistical Computing, 8, + ! 305. + + ! The only difference between the specific functions is the rank + ! of the first argument and the equal rank of the result. + + ! real, intent(in), rank >= 1:: f ! (n, ...) cell averages, n must be >= 1 + ! real, intent(in):: x(:) ! (n + 1) cell edges + ! real slopes, same shape as f ! (n, ...) + + module procedure slopes1, slopes2, slopes3, slopes4 + end interface + + private + public slopes + +contains + + pure function slopes1(f, x) + + real, intent(in):: f(:) + real, intent(in):: x(:) + real slopes1(size(f)) + + ! Local: + integer n, i + real xc(size(f)) ! (n) cell centers + real h + + !------------------------------------------------------ + + n = size(f) + forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. + slopes1(1) = 0. + slopes1(n) = 0. + + do i = 2, n - 1 + if (f(i) >= max(f(i - 1), f(i + 1)) .or. f(i) & + <= min(f(i - 1), f(i + 1))) then + ! Local extremum + slopes1(i) = 0. + else + ! (f(i - 1), f(i), f(i + 1)) strictly monotonous + + ! Second order slope: + slopes1(i) = (f(i + 1) - f(i - 1)) / (xc(i + 1) - xc(i - 1)) + + ! Slope limitation: + h = abs(x(i + 1) - xc(i)) + slopes1(i) = sign(min(abs(slopes1(i)), abs(f(i + 1) - f(i)) / h, & + abs(f(i) - f(i - 1)) / h), slopes1(i)) + end if + end do + + end function slopes1 + + !************************************************************* + + pure function slopes2(f, x) + + real, intent(in):: f(:, :) + real, intent(in):: x(:) + real slopes2(size(f, 1), size(f, 2)) + + ! Local: + integer n, i, j + real xc(size(f, 1)) ! (n) cell centers + real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) + + !------------------------------------------------------ + + n = size(f, 1) + forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. + + forall (i = 2:n - 1) + h(i) = abs(x(i + 1) - xc(i)) + delta_xc(i) = xc(i + 1) - xc(i - 1) + end forall + + do j = 1, size(f, 2) + slopes2(1, j) = 0. + + do i = 2, n - 1 + if (f(i, j) >= max(f(i - 1, j), f(i + 1, j)) .or. & + f(i, j) <= min(f(i - 1, j), f(i + 1, j))) then + ! Local extremum + slopes2(i, j) = 0. + else + ! (f(i - 1, j), f(i, j), f(i + 1, j)) + ! strictly monotonous + + ! Second order slope: + slopes2(i, j) = (f(i + 1, j) - f(i - 1, j)) / delta_xc(i) + + ! Slope limitation: + slopes2(i, j) = sign(min(abs(slopes2(i, j)), & + abs(f(i + 1, j) - f(i, j)) / h(i), & + abs(f(i, j) - f(i - 1, j)) / h(i)), slopes2(i, j)) + end if + end do + + slopes2(n, j) = 0. + end do + + end function slopes2 + + !************************************************************* + + pure function slopes3(f, x) + + real, intent(in):: f(:, :, :) + real, intent(in):: x(:) + real slopes3(size(f, 1), size(f, 2), size(f, 3)) + + ! Local: + integer n, i, j, k + real xc(size(f, 1)) ! (n) cell centers + real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) + + !------------------------------------------------------ + + n = size(f, 1) + forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. + + forall (i = 2:n - 1) + h(i) = abs(x(i + 1) - xc(i)) + delta_xc(i) = xc(i + 1) - xc(i - 1) + end forall + + do k = 1, size(f, 3) + do j = 1, size(f, 2) + slopes3(1, j, k) = 0. + + do i = 2, n - 1 + if (f(i, j, k) >= max(f(i - 1, j, k), f(i + 1, j, k)) .or. & + f(i, j, k) <= min(f(i - 1, j, k), f(i + 1, j, k))) then + ! Local extremum + slopes3(i, j, k) = 0. + else + ! (f(i - 1, j, k), f(i, j, k), f(i + 1, j, k)) + ! strictly monotonous + + ! Second order slope: + slopes3(i, j, k) = (f(i + 1, j, k) - f(i - 1, j, k)) & + / delta_xc(i) + + ! Slope limitation: + slopes3(i, j, k) = sign(min(abs(slopes3(i, j, k)), & + abs(f(i + 1, j, k) - f(i, j, k)) / h(i), & + abs(f(i, j, k) - f(i - 1, j, k)) / h(i)), slopes3(i, j, k)) + end if + end do + + slopes3(n, j, k) = 0. + end do + end do + + end function slopes3 + + !************************************************************* + + pure function slopes4(f, x) + + real, intent(in):: f(:, :, :, :) + real, intent(in):: x(:) + real slopes4(size(f, 1), size(f, 2), size(f, 3), size(f, 4)) + + ! Local: + integer n, i, j, k, l + real xc(size(f, 1)) ! (n) cell centers + real h(2:size(f, 1) - 1), delta_xc(2:size(f, 1) - 1) ! (2:n - 1) + + !------------------------------------------------------ + + n = size(f, 1) + forall (i = 1:n) xc(i) = (x(i) + x(i + 1)) / 2. + + forall (i = 2:n - 1) + h(i) = abs(x(i + 1) - xc(i)) + delta_xc(i) = xc(i + 1) - xc(i - 1) + end forall + + do l = 1, size(f, 4) + do k = 1, size(f, 3) + do j = 1, size(f, 2) + slopes4(1, j, k, l) = 0. + + do i = 2, n - 1 + if (f(i, j, k, l) >= max(f(i - 1, j, k, l), f(i + 1, j, k, l)) & + .or. f(i, j, k, l) & + <= min(f(i - 1, j, k, l), f(i + 1, j, k, l))) then + ! Local extremum + slopes4(i, j, k, l) = 0. + else + ! (f(i - 1, j, k, l), f(i, j, k, l), f(i + 1, j, k, l)) + ! strictly monotonous + + ! Second order slope: + slopes4(i, j, k, l) = (f(i + 1, j, k, l) & + - f(i - 1, j, k, l)) / delta_xc(i) + + ! Slope limitation: + slopes4(i, j, k, l) = sign(min(abs(slopes4(i, j, k, l)), & + abs(f(i + 1, j, k, l) - f(i, j, k, l)) / h(i), & + abs(f(i, j, k, l) - f(i - 1, j, k, l)) / h(i)), & + slopes4(i, j, k, l)) + end if + end do + + slopes4(n, j, k, l) = 0. + end do + end do + end do + + end function slopes4 + +end module slopes_m