Index: LMDZ5/branches/testing/libf/misc/regr1_conserv_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/misc/regr1_conserv_m.F90 (revision 2787) +++ (revision ) @@ -1,358 +1,0 @@ -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: LMDZ5/branches/testing/libf/misc/regr1_lint_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/misc/regr1_lint_m.F90 (revision 2787) +++ (revision ) @@ -1,98 +1,0 @@ -! $Id$ -module regr1_lint_m - - ! Author: Lionel GUEZ - - implicit none - - interface regr1_lint - ! Each procedure regrids by linear interpolation. - ! The regridding operation is done on the first dimension of the - ! input array. - ! The difference betwwen the procedures is the rank of the first argument. - module procedure regr11_lint, regr12_lint - end interface - - private - public regr1_lint - -contains - - function regr11_lint(vs, xs, xt) result(vt) - - ! "vs" has rank 1. - - use assert_eq_m, only: assert_eq - use interpolation, only: hunt !!, polint - - real, intent(in):: vs(:) - ! (values of the function at source points "xs") - - real, intent(in):: xs(:) - ! (abscissas of points in source grid, in strictly monotonic order) - - real, intent(in):: xt(:) - ! (abscissas of points in target grid) - - real vt(size(xt)) ! values of the function on the target grid - - ! Variables local to the procedure: - integer is, it, ns - integer is_b ! "is" bound between 1 and "ns - 1" - - !-------------------------------------- - - ns = assert_eq(size(vs), size(xs), "regr11_lint ns") - - is = -1 ! go immediately to bisection on first call to "hunt" - - do it = 1, size(xt) - call hunt(xs, xt(it), is) - is_b = min(max(is, 1), ns - 1) -!! call polint(xs(is_b:is_b+1), vs(is_b:is_b+1), xt(it), vt(it)) - vt(it) = ((xs(is_b+1) - xt(it)) * vs(is_b) & - + (xt(it) - xs(is_b)) * vs(is_b+1)) / (xs(is_b+1) - xs(is_b)) - end do - - end function regr11_lint - - !********************************************************* - - function regr12_lint(vs, xs, xt) result(vt) - - ! "vs" has rank 2. - - use assert_eq_m, only: assert_eq - use interpolation, only: hunt - - real, intent(in):: vs(:, :) - ! (values of the function at source points "xs") - - real, intent(in):: xs(:) - ! (abscissas of points in source grid, in strictly monotonic order) - - real, intent(in):: xt(:) - ! (abscissas of points in target grid) - - real vt(size(xt), size(vs, 2)) ! values of the function on the target grid - - ! Variables local to the procedure: - integer is, it, ns - integer is_b ! "is" bound between 1 and "ns - 1" - - !-------------------------------------- - - ns = assert_eq(size(vs, 1), size(xs), "regr12_lint ns") - - is = -1 ! go immediately to bisection on first call to "hunt" - - do it = 1, size(xt) - call hunt(xs, xt(it), is) - is_b = min(max(is, 1), ns - 1) - vt(it, :) = ((xs(is_b+1) - xt(it)) * vs(is_b, :) & - + (xt(it) - xs(is_b)) * vs(is_b+1, :)) / (xs(is_b+1) - xs(is_b)) - end do - - end function regr12_lint - -end module regr1_lint_m Index: LMDZ5/branches/testing/libf/misc/regr3_lint_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/misc/regr3_lint_m.F90 (revision 2787) +++ (revision ) @@ -1,100 +1,0 @@ -! $Id$ -module regr3_lint_m - - ! Author: Lionel GUEZ - - implicit none - - interface regr3_lint - ! Each procedure regrids by linear interpolation. - ! The regridding operation is done on the third dimension of the - ! input array. - ! The difference betwwen the procedures is the rank of the first argument. - module procedure regr33_lint, regr34_lint - end interface - - private - public regr3_lint - -contains - - function regr33_lint(vs, xs, xt) result(vt) - - ! "vs" has rank 3. - - use assert_eq_m, only: assert_eq - use interpolation, only: hunt - - real, intent(in):: vs(:, :, :) - ! (values of the function at source points "xs") - - real, intent(in):: xs(:) - ! (abscissas of points in source grid, in strictly monotonic order) - - real, intent(in):: xt(:) - ! (abscissas of points in target grid) - - real vt(size(vs, 1), size(vs, 2), size(xt)) - ! (values of the function on the target grid) - - ! Variables local to the procedure: - integer is, it, ns - integer is_b ! "is" bound between 1 and "ns - 1" - - !-------------------------------------- - - ns = assert_eq(size(vs, 3), size(xs), "regr33_lint ns") - - is = -1 ! go immediately to bisection on first call to "hunt" - - do it = 1, size(xt) - call hunt(xs, xt(it), is) - is_b = min(max(is, 1), ns - 1) - vt(:, :, it) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b) & - + (xt(it) - xs(is_b)) * vs(:, :, is_b+1)) / (xs(is_b+1) - xs(is_b)) - end do - - end function regr33_lint - - !********************************************************* - - function regr34_lint(vs, xs, xt) result(vt) - - ! "vs" has rank 4. - - use assert_eq_m, only: assert_eq - use interpolation, only: hunt - - real, intent(in):: vs(:, :, :, :) - ! (values of the function at source points "xs") - - real, intent(in):: xs(:) - ! (abscissas of points in source grid, in strictly monotonic order) - - real, intent(in):: xt(:) - ! (abscissas of points in target grid) - - real vt(size(vs, 1), size(vs, 2), size(xt), size(vs, 4)) - ! (values of the function on the target grid) - - ! Variables local to the procedure: - integer is, it, ns - integer is_b ! "is" bound between 1 and "ns - 1" - - !-------------------------------------- - - ns = assert_eq(size(vs, 3), size(xs), "regr34_lint ns") - - is = -1 ! go immediately to bisection on first call to "hunt" - - do it = 1, size(xt) - call hunt(xs, xt(it), is) - is_b = min(max(is, 1), ns - 1) - vt(:, :, it, :) = ((xs(is_b+1) - xt(it)) * vs(:, :, is_b, :) & - + (xt(it) - xs(is_b)) * vs(:, :, is_b+1, :)) & - / (xs(is_b+1) - xs(is_b)) - end do - - end function regr34_lint - -end module regr3_lint_m Index: LMDZ5/branches/testing/libf/misc/regr_conserv_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/misc/regr_conserv_m.F90 (revision 2839) +++ LMDZ5/branches/testing/libf/misc/regr_conserv_m.F90 (revision 2839) @@ -0,0 +1,388 @@ +MODULE regr_conserv_m + + USE assert_eq_m, ONLY: assert_eq + USE assert_m, ONLY: assert + USE interpolation, ONLY: locate + + IMPLICIT NONE + +! Purpose: Each procedure regrids a piecewise linear function (not necessarily +! continuous) by averaging it. This is a conservative regridding. +! The regridding operation is done along dimension "ix" of the input +! field "vs". Input are positions of cell edges. +! Remarks: +! * The target grid should be included in the source grid: +! no extrapolation is allowed. +! * Field on target grid "vt" has same rank, slopes and averages as "vs". +! * If optional argument "slope" is not given, 0 is assumed for slopes. +! Then the regridding is first order instead of second. +! * "vs" and "vt" have the same dimensions except Nr. ix (ns for vs, nt for vt) + +! Argument Type Description +!------------------------------------------------------------------------------- +! INTEGER, INTENT(IN) :: ix Scalar dimension regridded <=rank(vs) +! REAL, INTENT(IN) :: vs(*) Rank>=1 averages in source grid cells +! REAL, INTENT(IN) :: xs(:) Vector(ns+1) edges of source grid, asc. order +! REAL, INTENT(IN) :: xt(:) Vector(nt+1) edges of target grid, asc. order +! REAL, INTENT(OUT) :: vt(*) Rank>=1 regridded field +! [REAL, INTENT(IN) :: slope] Rank>=1 slopes in source grid cells + + INTERFACE regr_conserv + ! The procedures differ only from the rank of the input/output fields. + MODULE PROCEDURE regr1_conserv, regr2_conserv, regr3_conserv, & + regr4_conserv, regr5_conserv + END INTERFACE + + PRIVATE + PUBLIC :: regr_conserv + +CONTAINS + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr1_conserv(ix, vs, xs, xt, vt, slope) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:) + REAL, INTENT(IN) :: xs(:), xt(:) + REAL, INTENT(OUT) :: vt(:) + REAL, OPTIONAL, INTENT(IN) :: slope(:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt + REAL :: co, idt + LOGICAL :: lslope +!------------------------------------------------------------------------------- + lslope=PRESENT(slope) + CALL check_size(ix,SHAPE(vs),SHAPE(vt),xs,xt,ns,nt) + is=locate(xs,xt(1)) !--- 1<= is <= ns (no extrapolation) + vt(:)=0. + DO it=1, nt; idt=1./(xt(it+1)-xt(it)) + IF(xt(it+1)<=xs(is+1)) THEN !--- xs(is)<=xt(it)<=xt(it+1)<=xs(is+1) + CALL mean_lin(xt(it),xt(it+1)) + ELSE + CALL mean_lin(xt(it),xs(is+1)); is=is+1 + DO WHILE(xs(is+1)=1.AND.ix<=rk,msg) + ii=0 + DO is=1,rk; IF(is==ix) CYCLE + WRITE(msg,'(a,i0)')TRIM(sub)//" n",is + n=assert_eq(svs(is),svt(is),msg) + ii=ii+1 + END DO + ns=assert_eq(svs(ix),SIZE(xs)-1,TRIM(sub)//" ns") + nt=assert_eq(svt(ix),SIZE(xt)-1,TRIM(sub)//" nt") + + !--- Quick check on sort order: + CALL assert(xs(1)=1 source grid field values +! REAL, INTENT(IN) :: xs(:) Vector(ns) centers of source grid, asc. order +! REAL, INTENT(IN) :: xt(:) Vector(nt) centers of target grid, asc. order +! REAL, INTENT(OUT) :: vt(*) Rank>=1 regridded field + + INTERFACE regr_lint + ! The procedures differ only from the rank of the input/output fields. + MODULE PROCEDURE regr1_lint, regr2_lint, regr3_lint, & + regr4_lint, regr5_lint + END INTERFACE + + PRIVATE + PUBLIC :: regr_lint + +CONTAINS + + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr1_lint(ix, vs, xs, xt, vt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:) + REAL, INTENT(IN) :: xs(:) + REAL, INTENT(IN) :: xt(:) + REAL, INTENT(OUT) :: vt(:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1" + REAL :: r +!------------------------------------------------------------------------------- + CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt) + is = -1 ! go immediately to bisection on first call to "hunt" + DO it=1,SIZE(xt) + CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1) + r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb)) + vt(it)=r*vs(isb)+(1.-r)*vs(isb+1) + END DO + +END SUBROUTINE regr1_lint +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr2_lint(ix, vs, xs, xt, vt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:,:) + REAL, INTENT(IN) :: xs(:) + REAL, INTENT(IN) :: xt(:) + REAL, INTENT(OUT) :: vt(:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1" + REAL :: r +!------------------------------------------------------------------------------- + CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt) + is = -1 ! go immediately to bisection on first call to "hunt" + DO it=1,SIZE(xt) + CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1) + r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb)) + IF(ix==1) vt(it,:)=r*vs(isb,:)+(1.-r)*vs(isb+1,:) + IF(ix==2) vt(:,it)=r*vs(:,isb)+(1.-r)*vs(:,isb+1) + END DO + +END SUBROUTINE regr2_lint +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr3_lint(ix, vs, xs, xt, vt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:,:,:) + REAL, INTENT(IN) :: xs(:) + REAL, INTENT(IN) :: xt(:) + REAL, INTENT(OUT) :: vt(:,:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1" + REAL :: r +!------------------------------------------------------------------------------- + CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt) + is = -1 ! go immediately to bisection on first call to "hunt" + DO it=1,SIZE(xt) + CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1) + r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb)) + IF(ix==1) vt(it,:,:)=r*vs(isb,:,:)+(1.-r)*vs(isb+1,:,:) + IF(ix==2) vt(:,it,:)=r*vs(:,isb,:)+(1.-r)*vs(:,isb+1,:) + IF(ix==3) vt(:,:,it)=r*vs(:,:,isb)+(1.-r)*vs(:,:,isb+1) + END DO + +END SUBROUTINE regr3_lint +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr4_lint(ix, vs, xs, xt, vt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:,:,:,:) + REAL, INTENT(IN) :: xs(:) + REAL, INTENT(IN) :: xt(:) + REAL, INTENT(OUT) :: vt(:,:,:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1" + REAL :: r +!------------------------------------------------------------------------------- + CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt) + is = -1 ! go immediately to bisection on first call to "hunt" + DO it=1,SIZE(xt) + CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1) + r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb)) + IF(ix==1) vt(it,:,:,:)=r*vs(isb,:,:,:)+(1.-r)*vs(isb+1,:,:,:) + IF(ix==2) vt(:,it,:,:)=r*vs(:,isb,:,:)+(1.-r)*vs(:,isb+1,:,:) + IF(ix==3) vt(:,:,it,:)=r*vs(:,:,isb,:)+(1.-r)*vs(:,:,isb+1,:) + IF(ix==4) vt(:,:,:,it)=r*vs(:,:,:,isb)+(1.-r)*vs(:,:,:,isb+1) + END DO + +END SUBROUTINE regr4_lint +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE regr5_lint(ix, vs, xs, xt, vt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: vs(:,:,:,:,:) + REAL, INTENT(IN) :: xs(:) + REAL, INTENT(IN) :: xt(:) + REAL, INTENT(OUT) :: vt(:,:,:,:,:) +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: is, it, ns, nt, isb ! "is" bound between 1 and "ns - 1" + REAL :: r +!------------------------------------------------------------------------------- + CALL check_size(ix,SHAPE(vs),SHAPE(vt),SIZE(xs),SIZE(xt),ns,nt) + is = -1 ! go immediately to bisection on first call to "hunt" + DO it=1,SIZE(xt) + CALL hunt(xs,xt(it),is); isb=MIN(MAX(is,1),ns-1) + r=(xs(isb+1)-xt(it))/(xs(isb+1)-xs(isb)) + IF(ix==1) vt(it,:,:,:,:)=r*vs(isb,:,:,:,:)+(1.-r)*vs(isb+1,:,:,:,:) + IF(ix==2) vt(:,it,:,:,:)=r*vs(:,isb,:,:,:)+(1.-r)*vs(:,isb+1,:,:,:) + IF(ix==3) vt(:,:,it,:,:)=r*vs(:,:,isb,:,:)+(1.-r)*vs(:,:,isb+1,:,:) + IF(ix==4) vt(:,:,:,it,:)=r*vs(:,:,:,isb,:)+(1.-r)*vs(:,:,:,isb+1,:) + IF(ix==5) vt(:,:,:,:,it)=r*vs(:,:,:,:,isb)+(1.-r)*vs(:,:,:,:,isb+1) + END DO + +END SUBROUTINE regr5_lint +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +SUBROUTINE check_size(ix,svs,svt,nxs,nxt,ns,nt) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix, svs(:), svt(:), nxs, nxt + INTEGER, INTENT(OUT) :: ns, nt +!------------------------------------------------------------------------------- +! Local variables: + INTEGER :: rk, is + CHARACTER(LEN=80) :: sub, msg +!------------------------------------------------------------------------------- + rk=SIZE(svs) + WRITE(sub,'(a,2i0,a)')"regr",rk,ix,"_lint" + CALL assert(ix>=1.AND.ix<=rk,TRIM(sub)//": ix exceeds fields rank") + DO is=1,rk; IF(is==ix) CYCLE + WRITE(msg,'(a,i1)')TRIM(sub)//" n",is + CALL assert(svs(is)==svt(is),msg) + END DO + ns=assert_eq(svs(ix),nxs,TRIM(sub)//" ns") + nt=assert_eq(svt(ix),nxt,TRIM(sub)//" nt") + +END SUBROUTINE check_size +! +!------------------------------------------------------------------------------- + +END MODULE regr_lint_m +! +!------------------------------------------------------------------------------- + Index: LMDZ5/branches/testing/libf/misc/slopes_m.F90 =================================================================== --- LMDZ5/branches/testing/libf/misc/slopes_m.F90 (revision 2787) +++ LMDZ5/branches/testing/libf/misc/slopes_m.F90 (revision 2839) @@ -1,226 +1,283 @@ -module slopes_m +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 + ! Extension / factorisation: David CUGNET + + IMPLICIT NONE + + ! Those 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. + + ! slope(ix,...) acts on ix th dimension. + + ! 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, ...) + INTERFACE slopes + MODULE procedure slopes1, slopes2, slopes3, slopes4, slopes5 + END INTERFACE + + PRIVATE + PUBLIC :: slopes + +CONTAINS + +!------------------------------------------------------------------------------- +! +PURE FUNCTION slopes1(ix, f, x) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: f(:) + REAL, INTENT(IN) :: x(:) + REAL :: slopes1(SIZE(f,1)) +!------------------------------------------------------------------------------- +! Local: + INTEGER :: n, i, j, sta(2), sto(2) + 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) + REAL :: fm, ff, fp, dx +!------------------------------------------------------------------------------- + n=SIZE(f,ix) + 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 + slopes1(:)=0. + DO i=2,n-1 + ff=f(i); fm=f(i-1); fp=f(i+1) + IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN + slopes1(i)=0.; CYCLE !--- Local extremum + !--- 2nd order slope ; (fm, ff, fp) strictly monotonous + slopes1(i)=(fp-fm)/delta_xc(i) + !--- Slope limitation + slopes1(i) = SIGN(MIN(ABS(slopes1(i)), & + ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes1(i) ) + END IF + END DO + +END FUNCTION slopes1 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +PURE FUNCTION slopes2(ix, f, x) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: f(:, :) + REAL, INTENT(IN) :: x(:) + REAL :: slopes2(SIZE(f,1),SIZE(f,2)) +!------------------------------------------------------------------------------- +! Local: + INTEGER :: n, i, j, sta(2), sto(2) + REAL, ALLOCATABLE :: xc(:) ! (n) cell centers + REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) + REAL :: fm, ff, fp, dx +!------------------------------------------------------------------------------- + n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-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 + slopes2(:,:)=0. + sta=[1,1]; sta(ix)=2 + sto=SHAPE(f); sto(ix)=n-1 + DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) + DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) + ff=f(i,j) + SELECT CASE(ix) + CASE(1); fm=f(i-1,j); fp=f(i+1,j) + CASE(2); fm=f(i,j-1); fp=f(i,j+1) + END SELECT + IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN + slopes2(i,j)=0.; CYCLE !--- Local extremum + !--- 2nd order slope ; (fm, ff, fp) strictly monotonous + slopes2(i,j)=(fp-fm)/dx + !--- Slope limitation + slopes2(i,j) = SIGN(MIN(ABS(slopes2(i,j)), & + ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes2(i,j) ) + END IF + END DO + END DO + DEALLOCATE(xc,h,delta_xc) + +END FUNCTION slopes2 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +PURE FUNCTION slopes3(ix, f, x) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + 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, sta(3), sto(3) + REAL, ALLOCATABLE :: xc(:) ! (n) cell centers + REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) + REAL :: fm, ff, fp, dx +!------------------------------------------------------------------------------- + n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-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 + slopes3(:,:,:)=0. + sta=[1,1,1]; sta(ix)=2 + sto=SHAPE(f); sto(ix)=n-1 + DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) + DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) + DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) + ff=f(i,j,k) + SELECT CASE(ix) + CASE(1); fm=f(i-1,j,k); fp=f(i+1,j,k) + CASE(2); fm=f(i,j-1,k); fp=f(i,j+1,k) + CASE(3); fm=f(i,j,k-1); fp=f(i,j,k+1) + END SELECT + IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN + slopes3(i,j,k)=0.; CYCLE !--- Local extremum + !--- 2nd order slope ; (fm, ff, fp) strictly monotonous + slopes3(i,j,k)=(fp-fm)/dx + !--- Slope limitation + slopes3(i,j,k) = SIGN(MIN(ABS(slopes3(i,j,k)), & + ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes3(i,j,k) ) + END IF + END DO + END DO + END DO + DEALLOCATE(xc,h,delta_xc) + +END FUNCTION slopes3 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +PURE FUNCTION slopes4(ix, f, x) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + 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, m, sta(4), sto(4) + REAL, ALLOCATABLE :: xc(:) ! (n) cell centers + REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) + REAL :: fm, ff, fp, dx +!------------------------------------------------------------------------------- + n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-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 + slopes4(:,:,:,:)=0. + sta=[1,1,1,1]; sta(ix)=2 + sto=SHAPE(f); sto(ix)=n-1 + DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) + DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) + DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) + DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) + ff=f(i,j,k,l) + SELECT CASE(ix) + CASE(1); fm=f(i-1,j,k,l); fp=f(i+1,j,k,l) + CASE(2); fm=f(i,j-1,k,l); fp=f(i,j+1,k,l) + CASE(3); fm=f(i,j,k-1,l); fp=f(i,j,k+1,l) + CASE(4); fm=f(i,j,k,l-1); fp=f(i,j,k,l+1) + END SELECT + IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN + slopes4(i,j,k,l)=0.; CYCLE !--- Local extremum + !--- 2nd order slope ; (fm, ff, fp) strictly monotonous + slopes4(i,j,k,l)=(fp-fm)/dx + !--- Slope limitation + slopes4(i,j,k,l) = SIGN(MIN(ABS(slopes4(i,j,k,l)), & + ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes4(i,j,k,l) ) + END IF + END DO + END DO + END DO + END DO + DEALLOCATE(xc,h,delta_xc) + +END FUNCTION slopes4 +! +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +! +PURE FUNCTION slopes5(ix, f, x) +! +!------------------------------------------------------------------------------- +! Arguments: + INTEGER, INTENT(IN) :: ix + REAL, INTENT(IN) :: f(:, :, :, :, :) + REAL, INTENT(IN) :: x(:) + REAL :: slopes5(SIZE(f,1),SIZE(f,2),SIZE(f,3),SIZE(f,4),SIZE(f,5)) +!------------------------------------------------------------------------------- +! Local: + INTEGER :: n, i, j, k, l, m, sta(5), sto(5) + REAL, ALLOCATABLE :: xc(:) ! (n) cell centers + REAL, ALLOCATABLE :: h(:), delta_xc(:) ! (2:n-1) + REAL :: fm, ff, fp, dx +!------------------------------------------------------------------------------- + n=SIZE(f,ix); ALLOCATE(xc(n),h(2:n-1),delta_xc(2:n-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 + slopes5(:,:,:,:,:)=0. + sta=[1,1,1,1,1]; sta(ix)=2 + sto=SHAPE(f); sto(ix)=n-1 + DO m=sta(5),sto(5); IF(ix==5) dx=delta_xc(m) + DO l=sta(4),sto(4); IF(ix==4) dx=delta_xc(l) + DO k=sta(3),sto(3); IF(ix==3) dx=delta_xc(k) + DO j=sta(2),sto(2); IF(ix==2) dx=delta_xc(j) + DO i=sta(1),sto(1); IF(ix==1) dx=delta_xc(i) + ff=f(i,j,k,l,m) + SELECT CASE(ix) + CASE(1); fm=f(i-1,j,k,l,m); fp=f(i+1,j,k,l,m) + CASE(2); fm=f(i,j-1,k,l,m); fp=f(i,j+1,k,l,m) + CASE(3); fm=f(i,j,k-1,l,m); fp=f(i,j,k+1,l,m) + CASE(4); fm=f(i,j,k,l-1,m); fp=f(i,j,k,l+1,m) + CASE(5); fm=f(i,j,k,l,m-1); fp=f(i,j,k,l,m+1) + END SELECT + IF(ff>=MAX(fm,fp).OR.ff<=MIN(fm,fp)) THEN + slopes5(i,j,k,l,m)=0.; CYCLE !--- Local extremum + !--- 2nd order slope ; (fm, ff, fp) strictly monotonous + slopes5(i,j,k,l,m)=(fp-fm)/dx + !--- Slope limitation + slopes5(i,j,k,l,m) = SIGN(MIN(ABS(slopes5(i,j,k,l,m)), & + ABS(fp-ff)/h(i),ABS(ff-fm)/h(i)),slopes5(i,j,k,l,m) ) + END IF + END DO + END DO + END DO + END DO + END DO + DEALLOCATE(xc,h,delta_xc) + +END FUNCTION slopes5 +! +!------------------------------------------------------------------------------- + +END MODULE slopes_m